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Stability Derivatives for a Slender Missile 
with Application to a 
Wing-Bodv-Vertical-Tail Configuration 


ARTHUR E. BRYSON, JR." 


Hughes Research and Development Laboratories 


ABSTRACT 


Theoretical stability derivatives are determined for a slender 
missile through use of the concept of inertia coefficients of its 
@oss section. The determined by the 
Velocity potentials for two-dimensional motion of the cross section 


inertia coefficients are 


intwo mutually perpendicular directions and for rotation about 
an axis perpendicular to the cross section. The inertia coeffi 
tients are derived for the cross section of a configuration con 
sisting of a slender body of revolution with low aspect ratio wings 


and unequal vertical tails 


INTRODUCTION 


Mo: Jones,” AND SPREITER® have shown how to 
estimate the aerodynamic forces on slender fuse- 
lages, wings, and wing-fuselage combinations by con- 
Midering the flow in planes perpendicular to the wing or 
fuselage axis to be uninfluenced by the change in the 
component of fluid velocity parallel to this axis. Munk 
and Jones made use of the concept of the apparent ad- 
ditional mass of the cross section in a two-dimensional 
low to find the lift and moment due to angle of attack. 
Mhis concept is enlarged upon here to include all six 
of the inertia coefficients of the cross section, utilizing 
them to find all the stability derivatives (except the 
axial force derivatives) of the three-dimensional con- 
figuration. 

Consideration is limited here to configurations where 
the aft end of the fuselage and all wing trailing edges lie 
ima plane perpendicular to the long axis of the configu- 
Mation. Extension to cases where trailing edges occur 
@head of this base plane can be made by considering 
the effect of the vorticity shed from these edges as 

hown by Jones. 

Received November 18, 1952 

* Research Engineer, Guided Missile Laboratories. 


INERTIA COEFFICIENTS OF THE CROSS SECTION 


Consider a slender body? moving through an infinite 
fluid medium. By ‘‘slender’’ we shall mean that the 
maximum dimension of any cross section of the body is 
small compared to the length of the body and that the 
dimensions of the cross section change slowly with dis- 
tance along the length of the body. Let x,y,z be axes 
fixed relative to the body with origin at the center of 
gravity of the body and with the x-axis in the direction 
of the length of the body (see Fig. 1). We consider 
only those motions of the body where the instantaneous 
velocity vector of any point in the body makes a small 


angle with the body x-axis. Let £,n,¢ be axes such that 


+t We shall use the term ‘‘body,”’ henceforth, to mean a general 
aerodynamic shape that includes wings, fuselages, wing-fuselage 


combinations, etc 





nea 





€=0 PLANE, FIXED 
Im FWID 
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Coordinate systems used in the analysis 
fixed in body. 
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the € = 0 plane is a plane at rest with respect to the induced by a two-dimensional body (the three-dime - 
fluid far away from the body and such that the é-axis is sional body's cross section) whose velocities and Shanf a = 
parallel to the x-axis at the instant under consideration. change with time. 3 = 
Let the n,f-axes be the projections of the y,z-axes onto The kinetic energy of the fluid in the 7-¢ plane (q - 
the = 0 plane. tually a kinetic energy per unit distance normal to th _~ 
We approximate the three-dimensional fluid flow  7-¢ plane) is given by‘ 7 
around the body by considering the fluid in the 7-¢ z — 
plane to be incompressible and the fluid motion to be T’ = 5 pSlv, w, pD)[A] - 

two-dimensional—1.e., uninfluenced by any part of the is ms '* 
body or fluid which does not lie in the y-¢ plane. The - 4 
motion in the n-¢ plane is thus a nonstationary motion where the symmetric matrix [A ] is given by in 
alone. 

The 

shape 


ds, 


I r oO 3 O¢2 de l 7 Ov; .m fore, 
viet ae a" 
A 12 Aj; 
5 l Oye O¢2 l O¢; 
i4i =ite, An Zalo ~ 5 g #15, ds, 25, ds, D ae ds IE since 
A 13) A 23) A 33 
l 03 l 3 3 
D ls, s ~ dd; 
af Pay pos 


7) 
n 











and w = velocity of the cross section in the ¢ di 
. ; rection D 
S = a reference area . : 
p = angular velocity of the cross section al?S 
D = a reference length about the £-axis 
p = fluid density f ( )ds = the line integral around the boundary oi 
¢1 = potential due to a unit velocity of the the cross section 
cross section in the 7 direction The dimensionless quantities A,, we shall call the 
2 = potential due to a unit velocity of the inertia coefficients of the cross section. 


cross section in the ¢ direction 
THE ForcES ACTING ON THE Bopy Cross SECTION 


¢3 = potential due to a unit angular velocity 
of the cross section about the -axis [he aerodynamic sideforce, downforce, and rolling 
v = velocity of the cross section in the 7 di- moment per unit length on the body at any given cross 
rection section are given by* 


mf ie) fH 


Ov Ow 
z) ad j-er” “ or” 4 9 0 
A — _ 2 a WwW —_ Vv { 
dt| Ow t Ov 



































or” 4 or’ o7” 
| ry | ae | Dew | 


The yawing moment and pitching moment per unit (2) by a change in the shape of the cross section. 


length are given by Furthermore, the motion of the cross section can change 
n! x’ in two ways: | |) by a change in the linear velocity of 
| RA = wa (4) the e.g. of the body and (2) by a change in the angular qT 
velocity of the body. Thus, 
Now 7” is a function of v, w, p and the six inertia coef | By 
ficients. Thus 7” can change with time in two ways: v= BU + rxl (5) | over 


(1) by a change in the motion of the cross section and w= al’ — qxf 
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a? a a, = 
where o7 To [ do 
qa = angle of attack of the body Lin or ~ es 
ae ne ‘linear velocities aan is 
3 = angle of sideslip of the body f , aie d 01 — = | + pS[A] dw (7) 
_ . uw * ra ‘ 
q = pitching angular velocity of di Ow dt dt 
he body si 1 o7” d 
‘ git : .)angular velocities pD D P 
y = yawing angular velocity °f| | D op | — dt | 
the body . : . : 
. From Eqs. (5) and (6) we see that 
U = forward velocity of the body z a a - 
dv — ; ; , 
\]l five of these quantities are functions of the time dt Pe ee ee ae 
alone. lu 
+ . - . - . w . T ry * 7 
The inertia coefficients are functions of x, since the dl aU + al — gx+ Uq (8) 
shape of the cross section changes with x. There- ; 
; ap 
fore, D Dp 
dt | a f | 
(d/dt) = (0/dt) — U(0/Ox) (6) ie 
Substituting Eqs. (5), (6), (7), and (8) into Eq. (3), 
Since [A] is a symmetric matrix, we have, finally, 
7 [ BD UD *D? x | rD x 
ly - + B— —  s + — 
l U? i? Dp U D 
o 
D ; aD UD gD? x + gD x 
—— Z’ = — [A] : a a : ; x [Al] a — © — 
pS l Ol U. D o D U D 
I’ pD* pD 
LD | U? J oe 
F % rD x | 
a 
( rD =) 7 gD x pD \ gD x 
B+ : [EA] + (« — : [F] + —[G a- 9) 
UD U pI! y ll up| S 
pD 
L U 4 
l 0, 0, 0 force Z, the rolling moment /, the pitching moment m, 
iE} = 0, 0, 0 and the yawing moment » can be found to second order 
| An, Loo, A> in terms of a, B, p,q, 7, &, B, p, 9, 7, and U. 
r 0 0 We shall define force and moment coefficients as 
y follows: 
[F]} = -| 0, 0, 0 — ee 
An, Aj», A); oe a 
; ” A [Cy, Cz, € ne te => 
_ = a ee a D *xn/D Tae ce oe x 
, Ars, As», Ass 3 / ef 2, a( (11) 
IG] Ay Ap, Aj; (1 2) pl ‘5 re w/D D D D I 
0, Q), 9 4 
where x, = distance from the c.g. to the nose of the 
body and x, = distance from the c.g. to the base plane 


THe TotraL Forces AND MOMENTS ON THE Bopy 


By integrating Eqs. (9) and (4) with respect to x 
over the length of the body, the sideforce V, the down 


of the body. 


Integrating Eqs. (9) and (4) and using the definitions 


of Eq. (11), we have, finally, 
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= i rD Xp | [BD UD] Trp?) 
Cc io 8 
v UD l + Ba rr 
. ‘ gD X»y ~xn/D ‘) aD UD *xn/D x * gD? 
Cs = —2[A c : —2 [A | i( a — 2 [A | d( ) 
’ et FS 2 | VNo)lv 7° ; Dp“ \pD U? 
D )»D? 
C, P . = 0 
i. a a l a a U* al L J 
f y and [F 
Carry 
Pe ens \(6 rD 4 rE) + ( gD 4 iF) + pD iG] ( (= Pe all - rD x 7 
af ) A — UD [2 | Q@ UV D Uv 7 ( a a 2 f - ) A V D 
pD 
Ul 
[ = =6rD 
U 
gD 4 ale pl ) x ( qD 
es | F] G 1 = ( 
(« vp) +e lel “ND u| 
0 
; C } _ WD J 8 where 
‘ Uy D 
 —— qD X»y *xn/D ( =) \K; 
=—Ce i = 2 [A] —?2 [/ - 
Ce D a + U D / ook 1 | d D Qa 
C pD pD Ik. 
a al a U al L Us . 
rap rD UD) [ *D* ec 
a oe - and ya 
l l U* U’ analys 
*xn/D x x\|a@D gD UD *xn/D a\? (x gD? 
2 ai(*)a(5) = eae —2 ai(*) a(%) — ~ — 
«b/D D D l l U? oF cual D D Uy? 
pD* m Let 
[ pD 7 a pD rD | 


B, 
U1 t aes 
*xn/D x x D *xn/D x \* x »D qD 
Zz a (Ja ( ) Ae] —72 ic (*) a(*) ae = 5) 
‘ D D l ecraN D D U ; 
: Jj , All of 
(°°) 0 other 
L } 


_ “a ¢ 5 
Fro 














where (~) indicates the value in the base plane and the quantity Cy has no physical significance. 


Example: Consider the flat delta wing of span 6 and chord c. Here, 





0 Q, 0 
0 el ) / bi 
’ ai _ (Xn _ so 
{ bs where J ‘ LAS “> 
lc = £4 + & . 
0, 0 = 
L 64 SD? bd 





For simplicity we let S = (7 1)b? and D = 5; then, 








D? 


D? 


(12 
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7] T 7 [ b\? 7 
0. vg, 0 0, V, 0 0, -—-t<=], O 
b 
b\? ? ; 
[4] =]0, (;), 0 , [E) = |0, 0, 0| [G] = |0, 0, 0 (15) 
) 
0, O (55) O (5) QO 0, O ( 
~~ o y ‘igi 4 
and [F] vanishes. 
Carrying out the integrations indicated, we find, from Eqs. (12) and (13), 
2c pD bD gD 
a oe ins ie wee ee 
3D I Ud 
C. = ‘i , X gD 2c aD 2 -¢ UD 4K qD° 
a a Lee Sa oe 'U? 
1 pD 1 c pD? So < rD D DrD (16 
GC = - = ‘etal _* <2 : af — K, (a - 6%) nid = : : ‘ 
S$ 2 40 D U? 3 D l l os 


, Xp» l “) (*)' gD , (F ~~) _ UD _ gD? 
= — 7 _ a—2 -— kK -_— K — Kz | 
¢ 3 D) * p/ 1 ae uv) + Tp -U 


C, = Ki(pD/U)a — Kx(pD/U)(qD/U) 


where 
P l c 2 ? c Xp 
iki - ) ( 7 
6 \D 3D D 


Ik ” (2) (Sg) +25 (4) 
15 \D 3\D/ DD’ 38D\D 


These expressions contain all the terms of reference 5 and a few more. It is interesting to note that the sideforce 
and yawing moment terms that arise due to leading-edge suction are taken care of very simply in this method of 
analysis. 
THE STABILITY DERIVATIVES 
Let us denote derivatives in the following manner: 


( bet det beef 0 MM NO eM 


Qo 9 ) (e) fe) o re) re) fe) e) fe) 


Oa’ 0’ »D\’ D\’ _ (rD\’ w\' _ (BD\’ »D2\’ sD2\' 7D2\' UD (17) 
#952) 0 (0) (2) 0(22) 0(2) (MY 012) 0( 2) 0( | 
‘Ug l l l l UP U? U? u? 
All of the derivatives are evaluated ata = a,8 =p=q=r=p=qG=%=U=4&4= 8 = 0 (although some 
other equilibrium condition could, of course, have been specified). 
From Eqs. (12), (13), and (17) we see that we can obtain 55 stability derivatives. They are given as follows: 
C\ B» Cy a? Cy a An, Ay, Ais tn/D 0, Q, A 29 x 
Czz, Ce, Cz, = —2 Ap, Ao, Abs + Pa / 0, Q, = Aw d (=) ( 18) 
Cig, Cas C,,, é 113, A, Ags veiled A 11 A 22» 2A 12s A 13 
es i , X ral _ F D 2 : x (*) , 
—~L&Zry Zq eae “112, “129 — 20 . ( ( (19 
@ he ~ ) ) 
L ~—Cw, Cy ad A A m/D | A nu-Aw, 2A» f l 
[ abs ee Ce a i An, An, Ais Los "~/P TA, An, An q x\ 
lc -C,, —Cnp “DD {LAw, Aw, Ass ~ w/D LAw, Ax, Ass D 
*xn/D . ° 
0, 0, —Aw| * x 
» -- _ » 
a / i he 0, Pd D d (5) (20 
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Cis Cres Ci; ; 
Cy; —Cy, 
Cz: —Cz, ==—@ f 
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*xn/D ° = 

: A , oe x Xx 
12 9 ll» 12 ; 
“4 * 2 f 2 /D ra ‘ p@ (5) 21 
*y,/p | Au, Aw As x 

Ax, A», Ao d (*) (29 
*b/D A3, Aos, A33 
xn/D An, Ap x x 

Aj, Ao D d (=) (23 
adel Aj3, Ao 

./D 


Ani 2 *) 
, D . (= 


Cntr Cng 9 “/P TAn, Ar (=) i(*) 
, , ae 25 
Cm C mq — xb/D Ay, Ax D D 
Cy; *xn/D Ai. x 
Cz; = — 2a / Ao d (<) (96 
LCi ~#/? LAs 


= — 2a 


"xn 
/ xpb/D 


(27 


Aw] * i(*) 
Ax» | D D 


Thus, in general, for a body with a cross section with no symmetries, 14 integrals must be evaluated, 


Xn A ’ A 2» A 5 ° ia 0 . *y ) vo 
f j 0 : An. re d (5) / — Ay, Ap, Ais| * d (*) i" . = S) d (;) 
—xb/D 0. 0 a 1. D . xp/D 0, A», A»; D D ' —xp/D 0, Ax D D 
’ ’ £133 


For across section with one plane of symmetry (such as almost all present-day airplane and missile configurations), 


Ay = Ax» = 0, and, hence, only nine integrals must be evaluated. 
planes of symmetry, A,; also vanishes and only seven integrals must be evaluated. 


For a cross section with two perpendicular 
If the cross section is the same 


when rotated 90° (such as a cruciform missile configuration), then, also, Aj, = Ag and only four integrals must 


be evaluated. 


It is interesting that the only stability derivatives that require a knowledge of ¢3, the potential due to rolling, 


are C,, and C),. 


— Cz; 


Also, because of the symmetry of [A ], we find that 


Cie Cy Gy = Cy 
Cig Cys, Cma = Cz; 
= Cys; Cre Cz; 


Czzg oe on = ‘oe 9 
Cue = —Crar Cog = Cre 
Cz, ad -t C Yq c mB bai ( na 
and, in addition, when a = 0, Cig = Cy, and C,, = Cz,. 


TO A SLENDER WING-Bopy-VERTICAL- 
TAIL CONFIGURATION 


APPLICATION 


The method described above is applied here to find- 
ing the longitudinal and lateral stability derivatives of 
a configuration consisting of a slender body of revolu- 
tion with low aspect ratio wings and unequal low aspect 
ratio vertical tails (see Fig. 2). The cross section of the 
configuration is mapped conformally onto a circle in 
order to determine the velocity potentials for the low- 
speed two-dimensional motion of the cross section. The 
inertia coefficients of the cross section are then deter- 
mined, leading finally to the stability derivatives of the 
three-dimensional configuration. 


CONFORMAL MAPPING OF THE WING-Bopy-VERTICAL- 
TAIL Cross SECTION ONTO A CIRCLET 
It 


consists of a circle with two equal opposing fins and two 


The cross section of the body is shown in Fig. 5. 


unequal opposing fins at right angles to the equal fins 
The plane of the cross section will be called the w plane 
The conformal mapping of the exterior of the finned 
circle onto a circle is most easily seen in terms of three 
consecutive simpler mappings: 

(1) First the w plane is mapped onto the & plane by 
a Joukowsky transformation that ‘‘flattens’’ the circle, 


+ The author is indebted to R. H. Edwards for pointing out 


this transformation to him 
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BTABILITY 


§&= w+ (a?/w) 


(2) Next the “‘slit’’ & plane is mapped onto the Z 
plane by the transformation 


2 = 2 


(3) Finally, the ‘‘flat plate’’ along the real axis in the 
Z plane is mapped onto a circle in the ¢ plane by a trans- 
lation and another Joukowsky transformation, 


h—f (h+f)? 
2 16¢ 


Combining these three transformations, the mapping 
of the finned circle onto a circle is 


g*\3 h—- (kh +f)? /? 
o+(w+) -|* ene al (28) 
w Z 16¢ 


On the boundary of the contour, 


t= |(h+f)/4]e” \ 


as. iv(’ fy h TF os a) — 


(' —J ar . = cos 0) -c- sa*| 


where the plus sign holds for the right half of the w 
plane and the minus sign holds for the left half of the w 


(29) 


plane. 


The angles in the ¢ plane are given by 


2h — Vc? + 4a?) 


cos & = 1 — 
h+f 
cos 6 = 1 — [2(h — c)/(h + f)] 
cos % = 1 — [2f/(h + f)] (30) 
cos 6, = 1 — [2(f — c)/(h +f)] 
fase ig os ee Vc? + 4a’) | 


k+f 


The expressions for c, h, and f in terms of a, s, 4, and 
fs are 


c = s{l — (a?/s?)] 
| ay a*\* 
ve O-5) te0 48) 
S~ i,° 
a\* a)" 
qe (1-2) +e +2) 


It is easily seen that the cross section of the configura- 
tion considered can be specialized to twelve particular 
cases; these are shown in Fig. 4. The mapping can be 
easily extended to the case of a fuselage with elliptical 
cross section by first mapping the ellipse onto a circle 
by a Joukowsky transformation and then proceeding 
as above. 
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A 29 
the case of no vertical fins suffices to find A». 
has already been done by Spreiter.* 
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Fic. 2. Slender wing-body-vertical-tail configuration. 
THE INERTIA COEFFICIENTS OF THE WING-Bopy- 


VERTICAL-TAIL Cross SECTION 


The expressions for the inertia coefficients of the cross 


section are given in Eq. (2). 


Because of the plane of symmetry of the cross section 


under consideration here, we see that 


(32) 


Aw = Ao3 = (0 


The vertical fins do not affect the inertia coefficient 
since they lie along stream lines of the flow. Hence, 
This 


The result is 


A (3:3) 


. (=) (: a’ 4 a ) 
FrL D s* si 


The expression for A, is worked out in Appendix 


(1). The result is 


1/k+f\? _. : 2a\? 
> ( D ) [G(0, 64) + G(6;, 83) | —_ (=) (34) 


2 [ _ (= r : cos “)'| - 
Tv é 

— Xr — 
)K(e) | (* - .’) + (7 ——F r’) - 


cos A + cos u 
r mM 


2 cos = cos 
») » 


where 


G(A, 4) = 


— [K(k) — Eo(k)] X 
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Fic. 3. Conformal mapping of cross section onto a circle. plane ; 
and "~ leo 
Fig, k) = / = = incomplete ellip ( 
X . 0 V1 —k?*sin’ ¢ 
k = tan ; tan 9 (noteOQ <k < 1,sinceO SAS tic integral of first kind 
™OSuSar—d) K(k) = f # e) = complete elliptic integral ol | 
k’ = V1 —-FR first kind | This | 
short- 


E P “¢ b? oi ; A graph of the function G(A,u) is shown in Fig. 5, | 

1(¢, k) = V1— k* sin® ¢ (dg = incomplete - . : ae oT | the 

° " or incomplete and a table of values of the function is given in Table | he ge 
elliptic integral of second kind) A discussion of the function is given in Appendix (1). | damp 


The expressions for Aj3 and Az; cannot be given 1 


E(k) = B( ,) = complete elliptic integral of terms of tabulated functions. Phe expre ssion for A 
2 will not be given here, but Aj; is worked out in Appendix 
second kind (2) The result is 
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1/2a\-h+f.. \ ; 1 /h + f\3 . ; . ; 
( *) 2D (sin 0) + sin 6; — sin 0; — sin 64) — ( ‘) )2(sin 6) — sin 6, + sin 63; — sin 04) — 
ed Tv 


Av = ND 2D 
9 
= (sin* 0) — sin* 6, + sin*® 6; — sin*® 6,) — (cos 6) + cos 6; — cos 63; — cos 64) X 
0 
l ey : ’ , 
a (09 — 6; — 03 + O04) + (sin 26) — sin 26; — sin 203; + sin 204) | + 
T 
(cos 6) + cos 6; + cos 03 + cos 04)(cos 6) — cos 6; + cos 6; — cos 64) X 
Ss 
(ces 0) — cos 6; — cos 63 + cos 64) + 277(, 01, 83, 64) + 247(01, Oo, 44, as) (36) 
where 
H(a, 8, y,6) = H*(a, B, y, 6) + H*(6, y, 8, a) 


H*(a, 8, y, 6) = Fi V (cos 86 — cos a)(cos 6 — cos B)(cos 6 + cos y)(cos 6 + cos 6) cos 6 dé 
0 


The value of //* can easily be determined numerically for any particular cross section, since the integrand has no 


singularities. 
THE STABILITY DERIVATIVES FOR A WING-Bopy 2W 
VERTICAL-TAIL CONFIGURATION as a a’ Xp? 
goUme(1-=+-_\i1+— 
h. ao r° 


The stability derivatives for the equilibrium condi 


tona = a, 8 = p ae et eee B —~p=qd="' =< where W = weight of the missile, » = radius of gyration 
U = 0can be obtained by substituting the expressions oF the missile about the y-axis, and g = acceleration due 
i a Sa ete a eee a ; 

ior the inertia coefficients into Eqs. (13) through (27). to gravity. 


In general, the integrals involving Ay, and Aj; will 








have to be integrated numerically. The integrals in = eee tee 
volving A». can often be integrated in simple forms be- 
cause of the simplicity of the analytic expression for 


A». 


It is interesting that for a = 0 the only stability de nd 





rivatives that depend upon the distribution of the inertia 


coefficients with x are C Cag, Caps Curr Cue and all 








derivatives with respect to a, 8, p, g, 7, and U. All the 

other stability derivatives depend only on the values of (a) (b) (c) 
the inertia coefficients in the base plane. Furthermore, , 

the important quantities C,, — Cyg, and Cng + Cng 

depend only on the values of A; and Ag: in the base ee 


plane; in fact, (d) (e) (f) 


Car — Cag = — 2 (F) A, = (=) Cre 
7 D D ‘ 








| 





wa © 
Om + Cng = — 2 Ay = Cza 


D 
(9) (h) (i) 








This means, for example, that the damping of the 
short-period longitudinal oscillation depends only on 
the geometry of the base plane. If 7 is the time to 
damp to 1 /e times the undamped amplitude, then ———_ 


4W 


T = D? D\2 (j) (k) (1) 
gel’ | - ce ~ ( ) (Cng + Coa) 
| F Fic. 4. Special cases of general cross section 
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be Substi' 


Fic. 5. Plot of G(\,u). integr: 





Appendix (1) | ro 


INERTIA COEFFICIENT A}; 


Referring to Fig. 3, we wish to find the inertia coefficient (apparent additional mass coefficient) for motion of the 
cross section in the w, direction where w = w; + iw. The velocity potential for unit velocity in this direction is | 


¢1 = Rit —-1(¢ — w) + i[(h + f)?/16¢]} (1) | 


where Hence 
h’—f (A +f)??? P e**\' 
| — a o> —- =cC+\|(wt 
2 16¢ w 
On the boundary of the cross section ¢ = [(h + f)/4]e’’, so the potential there is 
¢: = [(A + f)/2] sin 0 — we where 
P h+f / F 9 
W, + WW. = + ; ~ [V (cos 8 — cos 4) (cos 6 + cos 6;) + V (cos 6 — cos 6;)(cos 6 + cos 63) } (2 


The inertia coefficient is given by 


] O¢) 4 9 
Ay = — ds = — e pidw (3 
. sf * on ™D* g ain where 
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TABLE |! 
The Function G(A,u) = G(u,d) 
\/u th 10° 20° 30° $0)° 50° 60° 70° 80° 90° 
0° 0) 
10 0.015 0.080 
20 0.059 0.074 0.117 
20) 0.129 0.148 0.185 0.250 
40) 0.220 0.234 0.273 0.334 0.413 
5) 0.325 0.338 0.374 0.431 0.508 0. 586 
60 0.437 0.449 0.482 0.534 0.599 0.673 0.750 
70 0.549 0.560 0.590 0.685 0.693 0.757 0.823 0.883 
g) 0.655 0.665 0.690 0.730 0.779 0.833 0.886 0.933 0.969 
90 0.750 0.757 0.779 0.812 0.853 0.897 0.937 0.970 0.992 1.0 
100 0.829 0.836 0.853 0.880 0.912 0.945 0.973 0.992 1.0 
110 0.891 0.897 0.910 0.931 0.955 0.977 0.993 1.0 
120 0.937 0.942 0.951 0.966 0.982 0.994 1.0 
130 0.968 0.971 0.978 0.987 0.995 1.0 
140 0.986 0.988 0.992 0.996 1.0 
150 0.995 0.996 0.998 1.0 
160 0.999 0.999 1.0 
170 0.999 1.0 
180 1.0 
From Eq. (2) we see that 
V (cos 6 — cos &)(cos 6 + cos 0;) + V (cos 6 — cos 6;)(cos 6 + cos 43); 0 < 6 < A% 
(4) 


Wy 


= + , V (cos 8 — cos 6;)(cos 6 + cos 03); 0% <O0< 0; r—-0; [0S wt — 


Also, it is easily seen that 


$ w. dw, = ma" (5) 


Substituting Eqs. (2) and (4) into Eq. (3), integrating by parts, and changing integration variables in the last two 


integrals, we have 


l/h + f\? Bo 
Ay = V (cos 8 — cos 0))(cos 86 + cos 64) cos 6 dé + 
7 D 0 


*@ P 63 
| V (cos 8 — cos 6;)(cos 6 + cos 63) cos 6 dé + f V (cos 8 — cos 63)(cos 6 + cos 6,;) cos 6 dé + 
0 0 


7A, a 2 
J V (cos 8 — cos 64)(cos 8 + cos 6) cos @ ao] —2 (<) (6) 
0 





| 
| 
| Hence, 


~] 
~— 


1 f(h+ ) (“) 
Ay = G(6o, 44) G(0,, 03)] — 2 ( 
11 »( D | 4 + 1 3 } D 


where 


G(A,u) = G*(A,u) + G*(p,A) ] 
: 2 A , (8) 
GAs) = V (cos 6 — cos A)(cos 6 + cos p) cos 6 0} 

7 JO 


where 0 < \ < randO<uw<7-.. 











The function G(A,u) can be expressed in terms of 
complete and incomplete elliptic integrals by a lengthy 
manipulation, leading to Eq. (9) in the text.+ 

Note that G(A,u) = G(u,d). If uw = 0, then k = 1, 
k’ = 0, and 


G(A,0) = 1 — cos! (A/2) (9) 
Ifu = 7 — X, then k = 1, k’ = O, and 
G(A, r — A) = 1 (10) 
If \ = yu, then 
G(A,A) = sin? A (11) 


It is easy to show that the contour lines of constant G 
have the following slopes (see Fig. 5): 


{ 0; «nz = 0| 
ay/dp = <—1l; A} = gpg (12) 
wo; X of 


t The author is indebted to H. A. Linstone for working out this 


expression. 


Two particular cases are easily worked out in elemen 
tary functions. 


4b: 
2t\? 2 
Ai —= ( ) (: - - 
D F 


which checks Spreiter’s result of Eq. (7). 


(a) 4h = fe = ¢, corresponding to Fip 


a 


(b) a = 0, tg = 0, corresponding to Fig. 4k: 
4 (i) (i , = 42° 1 +} “) 
é = -— 2 Z ( 
“ pb} \ ty? 1, \ ty? ‘ 
hy’. & ; 
Ay = as-~ — ©, corresponding to a flat plat: 
D Ss 


ty - § , 
Ay, = 23 as | — o, corresponding to a flat plate 
1 


in presence of an infinite end plate on one 


end 


Appendix (2) 


INERTIA COE 


The expression for Aj; is 


Ai = —(1/SD) f ¢1(0¢3/0n) ds 
8 "A "3 —6s 
Ais = Di (- J ; ¢, W, dw, + J s YW, dw» 
where 
¢ = [((A+ f)/2] sind — w, 


© 
| 


: 1f/h+f\? 
= } 9 [v (cos 6 — cos 6))(cos 6 + cc 


for0 < 0 < O@and rt — 6, <O0< =; 


FFICIENT A}; 


"r | l 
_ / QW aw,)| 
O=x—0, 


»s 03) + V (cos 6 — cos @;)(cos 6 + cos 


: l/h +f? / 
— = [v (cos 8 — cos 0)(cos 6 + cos 04) + V (cos 6; — cos 8)(cos 6 + cos 63) 


é { 


9 


for 6, 3 6 < = 6s. 


Substituting Eq. (2) into Eq. (1), integrating by parts, and manipulating, we finally arrive at Eq. (10) of the | 


text. 
For the special case of no body —i.e., a = 0 
P HELV in Oy + sin 6 bss 
3;=- \sin & sin 6, — ~ sin* 0) — 
- 2r \ D a : 3 
l 
(cos 04 — cos 6) (09 — O05) + 
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Pee 
' (sin 26) — sin 20) — cos 49 cos 44 (sin A) + 


Eq. (10) simplifies greatly, since 6) = 6; and 0; = 6, for this case 


l 
 sin® 04 + 


: | 
sin aT 5) 
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Calculation Method for Three-Dimensional 
Rotationally Svmmetrical Laminar Boundary 
Lavers with Arbitrary Free-Stream Velocity 
and Arbitrary Wall-Temperature Variation’ 


ROBERT M. DRAKE, JR.+ 
University of California al Berkeley 


ABSTRACT 


A calculation method is presented enabling a rapid prediction 
of the distribution of the local skin-friction and heat-transfer 
coeficients to be made for rotationally symmetrical bodies having 
arbitrarily specified free-stream velocity distributions and _ sur- 
face-temperature variations for constant-property, incompres 
sible, laminar boundary-layer flows. The method utilizes the 
Mangler transformation as a basis, thus enabling the results cal 
culated by Hartree and Levy for the two-dimensional case to be 
used directly in the determination of boundary-layer quantities 
for the three-dimensional case. 

Example results using the proposed method are presented in 
the form of the boundary-layer momentum and displacement 
thicknesses, the skin-friction coefficients, and the heat-transfer 
coefficients for a three-dimensional half-body with constant 
surface temperature. These results, utilizing the Hartree ve 
locity profiles, are compared with similar results obtained by 
Scholkemeier utilizing the Pohlhausen ?-4 velocity profiles. The 
skin-friction coefficient for stagnation-point flow as obtained by 
the approximate method is compared with the exact solution 
for this region given by Homann to show excellent agreement. 


NOMENCLATURE 
A, B, C = constants 


= thermal diffusivity, k/pcp (ft.?/sec. ) 
specific heat at constant pressure (B.t.u./Ib.-°F.) 


(Cy = skin-friction coefficient, to/!/2p lo 2, dimensionless 
D = maximum diameter of half-body (ft.) 
F(Z) = solution of momentum equation, Eq. (5) 
(Z = solution of energy equation, Eq. (28) 
= acceleration of gravity (ft./sec.*) 
h = heat-transfer coefficient (B.t.u./hr.-ft.2-°F.) 
k = thermal conductivity (B.t.u./hr.-ft.-°F-.) 
L = constant in Mangler transform (ft.) 
m = exponent defining free-stream velocity variation 
n = normal coordinate (ft. 
= heat flux (B.t.u./hr.-ft.?) 
r = thermal recovery factor, Tae — 71/(U1?/2g/Cp) 
r = radius of body (ft.) 
= tangential coordinate (ft.) 
S) = distance from stagnation point along body surface 


(ft.) 
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T = temperature in boundary layer (°F. ) 

T aw = adiabatic wall temperature, Tg, = 7) + 7 (Ui? 
2g) Cp)( °F.) 

= surface temperature of body (°F.) 


7; = state temperature at edge of boundary layer (°F.) 
1 = stagnation temperature at edge of boundary layer 

( °F.) 

U, = velocity of the free stream at edge of boundary 
layer (ft./sec.) 

Uy = velocity of the onstream flow (ft./sec.) 

u = tangential velocity component (ft./sec.) 

v = normal velocity component (ft./sec.) 

FA = dimensionless normal coordinate 

P Aug = dimensionless displacement thickness 

Z; = dimensionless momentum thickness 

a = included angle of cone, deg. 

B = parameter relating to the free-stream velocity 


variation, 2m/(m + 1) 
exponent defining temperature variation of surface 


/ 

6 = boundary-layer thickness (ft.) 

6* = boundary-layer displacement thickness (ft.) 

m = absolute viscosity (Ib. sec./ft.*) 

v = kinematic viscosity (ft.*/sec.) 

7 = stream function in the boundary layer (ft.*/sec.) 

6 = boundary-layer momentum thickness (ft.) 

0, = boundary-layer thermal momentum thickness 
(ft.) 

To = shear stress at the surface (Ib./ft.? 

p = density (Ib. sec.?/ft.*) 

o = Prandtl Number, ucp/k = v/a, dimensionless 


Nore: Barred symbols indicate two-dimensional values; 


unbarred symbols indicate three-dimensional values 


INTRODUCTION 


roves OF THE laminar boundary-layer char- 
acteristics on rotationally symmetrical, three 
dimensional bodies is significant to the design of high- 
speed aircraft and missiles. This rotationally sym 
metrical flow has been investigated by Millikan,' Fréss 
ling,” Scholkemeier,* Tomotika,‘ Whitehead and Ca 
netti,> and others, with the ultimate aim of specifying 
distribution of the local skin-friction and heat-transfer 
coefficients on the surface of any arbitrarily shaped 
body. The usefulness of any so devised specification 
lies in its utility, as well as the accuracy, of the results. 
In practical instances, a variation of the surface tem 
perature generally occurs which has a marked effect 
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Ou 4 Ou U dU, O*" 
u v = 1 v ; (moment 
Os on ds On? — 
oT oT o°T 4 v (x) 
u v =a x (energy ) (9 
Os On On? Cp \On 6) 5 
Us F [O(rou)/Os] + [O(rov)/On] = 0 (continuity) (3 
Zz . — 
with the general boundary conditions, 
n=0, u=v=0, T = 7,(S) 
LLLILLLL LLL LIM LL , : 
Fic. 1. Coordinate system for rotationally symmetrical body. ee OMe Fee 





upon the heat-transfer rate, and, as a consequence, 
solutions of the momentum and energy equations ac- 
counting for the effect of variable wall temperature, as 
well as variable free-stream velocity, are necessary. 


All boundary-layer solutions in which the variation 
of wall temperature can be accounted for relate to con- 
stant property, incompressible flows with a free-stream 
velocity distribution that is expressible as a power func- 
tion of the meridional distance from the forward stag- 
nation point (l; = C;S”). Exact solutions of the 
corresponding energy equations are possible for a wall- 
temperature distribution of the same type (7) = 
A;S”). Such solutions for the momentum and energy 
equations for two-dimensional flows have been com- 
puted by Hartree® and Levy.’ 


Since for the specification of an arbitrary free-stream 
velocity and surface temperature distribution there is 
correspondence with exact solutions only at the stag- 
nation point, approximate solutions are resorted to for 
the regions downstream of the stagnation point. Such 
approximate methods for two-dimensional flows have 
been discussed by Seban.* 


Since rotationally symmetrical three-dimensional 
flows, wherein the boundary-layer thickness is much 
smaller than the body radius, may be represented by 
momentum and energy equations that are identical 
with the two-dimensional case (the three-dimensional 
character of the flow having its specification only in the 
continuity equation), the two-dimensional flow results 
are directly applicable to these particular three-dimen- 
sional flows. It is the purpose of this paper to present 
an approximate method for the calculation of laminar 
boundary-layer characteristics for rotationally sym- 
metrical three-dimensional flows that will yield a spec- 
ification the local skin-friction and heat-transfer 
coefficients for an arbitrarily specified free-stream ve- 
locity and surface-temperature variation. 


of 


ANALYSIS 
In the instance wherein the boundary-layer thick- 
ness is small when compared with the body radius the 
applicable equations for constant property, incompres- 


sible flow, are (Fig. 1) 


Momentum Equations 


Exact Solutzon. 
satisfied by a stream function of the type, 


The continuity equation may be 


1 Oy 
u= v= 


ry On 


1 Oy 
Yo Os 


which statements when used in combination with the 
Mangler transformation,? where barred variables are 
two-dimensional quantities, 


. "S Cry? Cri 
Qo = / L? ds, R= 7 


result in a statement of the momentum equation of the 


form 
r= o'y aE oy (‘) Oy 
L) da dnods L)/ as on?” \L/ on 
ad; 
U 
di 


which yields, upon setting Y = (C/L)y, the following 
equation, which is identical in form with the correspond- 
ing two-dimensional case: 

oy 


on? 


oy oy 
On OnOsS 


dl 1 
ds 


3,7, 
O*y 


oy 
ao yar 
on* 


OS 
Eq. (1) may be further transformed by virtue of the 
Hartree variables, wherein 


U; _ CoS” 


z= m+ 1 Ui — oe . Ui 
2 vS V2-8 vs 


yp =V2-—BV UWS F(Z) 
to yield the familiar equation for the ‘“‘wedge flows’’ for 
which solutions are available.® 


F’'"(Z) + F(Z)F’'(Z) + B {1 — F°(Z)} =0 (5 


The relationship between the two- and three-dimen- 
sional cases may be shown by considering the three- 
dimensional body as a right circular cone in analogy 
with the two-dimensional wedge. The cone may be 
represented by 7) = S sin a and the velocity distribution 


r ‘ ‘ + . ‘ ‘el . 
as U; = C38”. Now relating S to S, one obtains 
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LAMINAR 


2 ms Cc? s 
f= : f rr dS = — fj S’ sin? adS = 
 £ Je L* Jo 


j C? sin? a Lc 
a 6S C65 
and for the velocity distribution, 


, JC? sin? al” 


U = CoSm si Co + S3m = C3S™ 
2 3 Cf 
Therefore, by comparison, 

C C jC? sin? a\”™ é m ‘ 

$= 9 * F > m= ) 

lL: 3 f 3 = 

Further, since m = B/(2 — B) and m = B/(2 — 8), 
the following relationships exist: 

B = 8/(3 — B); B = 38/(8 + 1) (7) 


Thus, by virtue of the above transformation, exact 
solutions are available for three-dimensional axisym- 
metrical flows wherein the free-stream velocity dis- 
tribution may be expressed as a power function of the 
distance from the stagnation point. 

Approximate Solution.—It is possible in the three- 
dimensional axisymmetrical case to utilize the exact 
solutions for the ‘‘cone flows’ to form approximate 
solutions for such bodies with arbitrary free-stream 
velocity distributions in a form that utilizes the compu- 
tations already made for the ‘‘wedge flows.” 

The momentum equation may be integrated across 
the boundary layer to yield the expression 


: thn es f / 
rei u* dn re u dn 
iy, f° i A a nse, = 
nds \' Jo — 0 9 is = ' ds 
(o) 
v (S) 
On/ » 


| dry 


Eq. (8) may be written as 


7) _. dé _ dl, 
= U;? + (20 + 6*) U, 
p ds ds 


6U;* (9) 


ry ds 
wherein 


- f “(1 — Man; w= f (1 = *) an 
o U, U; 0 U, 


Eq. (9) for the rotationally symmetrical three-dimen- 
sional case differs from the corresponding statement for 
the two-dimensional case only in the third term on the 
right-hand side. Using the Holstein-Bohlen form pa- 
rameter, (@2/v) (dU’,/dS), Eq. (9) may be written as 


™ 8 cig B+ (245) eet . 
uu, oy dS . 6 vy dS 
| Iro 6° : 
-—U,; (39) 
ro dS v 
Writing Eq. (10) in terms of the variable 6°/v gives 
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U, d(e? v) U; dro (*) _ 

2 dS ro AS \v 
{to 0 (2 “) dl i 
< — {2 —? 11 
lu U; as 


In the instance wherein the free-stream velocity dis- 
tribution may be expressed as a power function of the 
distance along the body LU; = C;S” and the dimension- 
less momentum thickness is defined as 


“ Fin . 

Z, = (6/V2 — B) VU;/vS 
the right-hand side of Eq. (11) is a function of the form 
parameter, which in this case is 

(@?/v) (dUi/dS) = Bz," 
By means of the Mangler transformations, the left- 

hand side of Eq. (11) yields the relationship 
1 U, de? 


U; dro ° 
(0°/v) = (12) 


U; d(6 v) 
ro dS 2 » d§ 


> ¢5 


Seban*’ has shown that this expression may be expressed 


as a linear function of the form parameter BZ. Thus, 
U, d(6?/v) U, dro, 1 U, de? 7 
2 dS nes” ' 06 Cc OS 

—ABpZ,?>+ B (13) 


where A and B are constants. It can be shown also 
that the form parameter is invariant in the trans- 
formation, so that 


BZ? = BZ? (14) 


Thus the constants determined for the two-dimensional 


case are applicable hereto.| Then one can write, from 


Eq. (13), 


+ B 


(15) 
v 


@? 6? dU 

d —A > 
v 4 ou") _ v dS 
dS ry dS _ U; 


which when rearranged to be in dimensionless form is 


6 _D . ro 
r (ry y ) ,2 4 (5) 
+ 


d(S/D) \n/D d (S/D) 


Ud Ui \ /0 : 
A = (= : Uy = B(=*) wo 
U; : f D v U; 
d(S/D) 
The solution of Eq. (16) yields the expression 
(? Pay 7 . ("*) 9 (2 ) ge D (“) 7 
D y D Ue 0 D 
U,\4-" (S 
(22) or 
U. D 


t The value of the constant A was changed slightly from the 


value given in reference 8. Change in final results was negligible. 
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which gives the momentum thickness as a function of 
the body shape and the free-stream Reynolds Number 


for regions away from the stagnation point. At the 
stagnation point S/D + 0; U,/U. > 0; m/D—>0; 
thus, 
0 aa) 
D v 7 
s 2 r\A-1 . 
ro \- l R 
f(s) ) 4G) 
lim 0 D U. D (18 
, o) 
S/D—-0 (") (; y 
D U. 


The limiting solution of Eq. (18) 1s 


VW B/(2 + A) 


(19) 
ia tS D) | 


The values taken for the constants are 


Jace 


rye 


0 < BZ;? < 0.09, 
—0.068 < BZ? < 0, A= 53 


Az §:2 
> 
) 


(19), the 
boundary-layer 


(17) and 


the 


Applying these constants to Eqs. 
final expressions are obtained for 


0.247 [yi 


momentum thickness. 


At stagnation point, 
5 ye d(U, ee) 
D - 


(20) 
1(.S/D) i 
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Functions relating to solutions of Eq. (5). 


< 0.09 


fordO<BLZ? 
fa) U sj ) 2 (; p ) 5.23 
: = ().44 
6 y : G U x 
7 (3) ( 1 \" i d (") 91) 
0 D l D 
ar 


and for —0.068 < BZ; <0 
y “ 2 ) 5.32 f D 
= 0.44 x 
(b) (Ce) 
(=) U 1 1.33 a 
(=) (>) @ 
D Be D 


lU.D 


ele 





(21), and (22), together with the knowledge 
and the body 
boundary-layer | 


Eqs. (20), 
of the free-stream velocity distribution 
the 

as a function of the distance S$ 


shape, yield a_ specification of 
momentum thickness 
The displacement thickness 6* 


along the body. ma\ 


be obtained by means of the relationship 
6*/0 = Z*/Z, = 2*/Z, 
where Z* is the dimensionless displacement thickness 
Z*/Z, = Z2*/Z 


is shown as a function of the form parameter on Fig. 2. | 
From this relationship and a knowledge of 6 the values 
may be computed. 


for 6* The form parameter is best 


calculated from the relationship 
UD)? d(U,/U.) 
vy ' d(S/D) 


j9 


BZ ,* = + 
\D 
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LAMINAR 


The values of the form parameter are necessary also for 
the selection of the values of the constants in the above 
equations. 


Skin-Friction Coefficient. The values of the skin- 
friction coefficients may be determined from the avail- 


able results. Thus, 


r — w(Ou/On)y 2 (=) (~) 
+ o(U.?/2) p(U.?/2) U.? L On] 


Since, 


l i= ( oe ( Sm 
u/U, = F(Z) 
as U,F'"(Z U — ? 
trib — = — : 1 = U;F'(Z) — 
on din VY2—8 vs 6 
Thus, 


: 2, _,, () ia 
C, = — F"’'(0) —* 
l 'g 6 
6 = (Cro i )0 
Therefore, 
GC F'(0) Z, 


~ vY oo ' ()rG D)VU.D/vy) “ 


i) 


The functions F’’(0) and Z, have been calculated and 
tabulated as functions of B, Fig. 2.6 '" The value of 6 
is determined from Eqs. (20), (21), and (22). 

Comparison with Homann Exact Solution for Stag- 
nation Flows.—-The results of the above analysis may 
be compared with the exact solution given by Homann!! 
for the skin-friction coefficient at a three-dimensional 
stagnation point. In this instance the free-stream ve- 
locity is a linear function of the distance S; thus m = 
lin the expression U7; = C3”. 

By Eqs. (6) and (7), therefore, m = '/3 and B = 1/9. 
For a value of 8B = '/. the values of F’’(0) = 0.9276 and 
Z, = 0.3502. The value of the momentum thickness 
at the stagnation point is given by Eq. (20), 


6 |UD 0.247 0.247 
me | y ge U i ¢ v.) 
y d(S/D) 4, S/D Jo 
Introducing these values into Eq. (23) gives 
G |U.D  F"O)zZ s (e =) : 
2 \ vy — 0.247 LU V S/D }y 
ee 
Loi} : Bi (24) 
l 


wherein $8; is the velocity gradient at the stagnation 
point defined by 


| 


U,/U. = Bi(S/D) 
The exact result given by Homann for symmetrical 
three-dimensional stagnation flow is 
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(C,/2) V U.D/v = 1.38172 (U,/U.) pi” (25) 


Comparison of Eqs. (24) and (25) shows excellent 


agreement. 


Energy Equation 


It will be assumed here that superposition of the 
solutions of the energy equations to account for the 
effects of mechanical energy dissipation holds on_a local 
basis as has been outlined in reference 8. The follow- 
ing development will consider first the case wherein 
dissipation is negligible. The necessary changes to 
account for dissipation will be indicated. 

Exact Solutions—The energy equation, Eq. (2), 
neglecting the dissipation term may be transformed 
by the Mangler transformation in the manner of the 
momentum equation to yield in the two-dimensional 
variables, Y, m, and §, 

wor oor oT 
- =a (26) 
On OS OS On On * 

To account for surface temperature variation, the 
temperature 7° is considered as a product, one part a 
function of S and the other a function of Z; there fol 


lows: 


g(O) = 1 


a 7T(.S)-¢ 7 * 
7 To(S)-2(Z); debae 


Transforming now to the Hartree variables as be- 


fore, one obtains 


wi 5 = § d1o “FF FF a. 
F'(Z)e(Z) (2 — B)AL= =) — F(Z)g"(Z) = - 2""(Z) 
To ds o 
(27) 
If 7) = AS’ = A;5S”, Eq. (27) reduces to the follow 


ing form: 


g’'(Z) + oF(Z)g'(Z) — oF(2 — B)F'(Z)g(Z) = 0 
(28) 


Eq. (28) has been solved for several values of the ex 
ponent 7 and for several values of the Prandtl Number, 
a, by Levy’ and for the special case of ¥ = 0 (i.e., con 
stant wall temperature) by Eckert.'® 

Approximate Solution.—-The energy equation may 
be integrated across the boundary layer to yield, for 
the rotationally symmetrical case, 


~ 


1 dry ae d en or 
ul’ dn + ul dn a 
ro ds Jo ds Jo on / >» 


(29) 


or when rearranged to be dimensionless in the inte 


grands, 
1 dry Ve TF d al " 
; .) dn + UT» x 
rm ds JO U; T% ds J0 
(" ) ~ ) 
— | dn a (30) 
U, 1% On/o 


UT» 
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If 


H -[ Pd 
= b 7 ae 


Eq. (30) becomes 


l dro 


UT + = aura « (o") 31) 
a awe OE 


Application of the Mangler transformation yields, for 
Eq. (31), the equivalent two-dimensional statement, 


d , oT 
UiT6,) = —- 32) 
d§ ( . " = ) ( 


Eq. (32) is satisfied on a local basis by the solutions of 
Eq. (28) within the range of the possible discrepancies 
previously noted regarding the variability of 8 as a 
parameter.’ It should be noted here, however, that 
the range of B needed for the three-dimensional calcu- 
lation is 0.5 < 8 < —0.20, whereas for the two-dimen- 
sional calculation the range of B needed is 1.0 < B < 
—0.20; therefore, whatever the error involved by neg- 
lecting dB /dS will be diminished in the three-dimensional 
case. 

Heat Transfer to the Boundary Layer.—The heat 
rate is given as follows:t 


oT : dZ (Cro 
=—k = —kA,S5%p'(0 ( ) 
= ) ew da\ L 
, Z, : Zi 
qo = — kA»S%g’(0) fs —kA;S"g’(0) ry (33) 


The momentum thickness is determined from Eqs. 
(20), (21), (22), while g’(0) and Z, are the two-dimen- 
sional functions corresponding to the applicable form 
parameter as determined from 6. g’(0) as is shown in 
Eq. (7) is relatively insensitive to 8 but depends most 
critically upon the value of 7. 


The heat rate determined from Eq. (33) may be ex- 
pressed in terms of the usual dimensionless parameters 
as follows: 

—g'(0)Z; 


- = (34) 
(0/D)V U..D/» 


(hD/k) 
VU.D/» 


It should be noted here that solutions for a power 
function variation in the wall heat rate may be com- 
puted on a similar basis to the method outlined above 
for the variation in the wall temperature. Briefly 


the method is as follows: From Eq. (33), 


+ The heat rate as calculated here does not satisfy the integral 
energy equation, Eq. (29), for the approximate case. It could 
be made to do so at the expense of not satisfying the momentum 
equation; however, the results of this approximation when ap- 
plied to experimental data, Eq. (12), is satisfactory, which fact 
when considered in conjunction with the rapidity of calculation 


makes the approximation acceptable. 


6 V1 iS 


0 = — kA3S%2’(0 ; 
q aid. oe 


(35 
wherein, by means of the transformation, it can be 
shown that 


6/0 = Cn/L = V3 VS/S (36) 
Inserting Eq. (36) into Eq. (35), one obtains, t 


S VU, vS 


qo = —kA3S"g'(0) V3 37 
, SV2—i (3 


From Eq. (37), 
qo = KS7vtla-) 2) a KS§7~(a-8) (2—8)] = KS? 


wherein p = y — [(1 — 8)/(2 — B)], the exponent for 
a power-function heat-rate distribution. Thus it can 
be seen that the problem may be presented on the basis 
of either heat-rate distribution or surface-temperature 
distribution. Knowledge of one allows a determination 
of the other. For instance, if the heat-rate distribution 
is given, the temperature distribution is defined by the 
exponents y, 6, or the corresponding two-dimensional 
values 7, 8, 


y= pt ([(1 — 8)/(2 — B)] 


As an example, if a constant heat rate is specified, p = 
0 and the temperature distribution exponent becomes 


y = (1 — B)/(2 — B) 


whereas if the temperature distribution is constant (iso- 
thermal wall), the heat rate must conform to values of 


p given by 
p = (8 — 1)/(2 — B) 


Other examples involving arbitrary distributions of the 
heat rate or surface temperature may be easily con- 
sidered. 

Effect of Frictional Heating—The assumption of 
superposition of solutions to account for the frictional 
heating effect on arbitrarily shaped bodies has been 
noted in reference 8 and to a great degree confirmed 
experimentally in reference 12. The method briefly is 
as follows: 

For wall-tem perature variation Ty) = Ty(S): The heat 
rate is computed from solutions of Eqs. (28) and (33 
(which neglect dissipation), with an effective temper- 


ature difference of 


t Eq. (37) states the analogy between flat plates and super- 
sonic cones (U,; = constant along surface) as originally stated by 
Handzsche and Wendt!*—i.e., that the local heat transfer on the 
surface of a supersonic cone was related to the local heat transfer 
on the surface of flat plate by the factor 4/3 when the distance 
from stagnation point to the point under consideration was the 
same foreach. Eq. (37) gives, for the condition S = Sand B = 9, 

) 


qo cone = V3 Jo flat plate 
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—— ia it U;? 
To —- Tew = To-—T1-17 = 
2gJcp 
- , uy . 
71> —7;,, + (1-7) (38) 
2g] cp 


wherein 
ture and r is the thermal recovery factor.'* 


7, is the free-stream stagnation tempera- 
The 
thermal recovery factor may be taken as the Vo for 
the region (1/2) < B < 0 as a satisfactory approxima- 
tion for the laminar boundary layer. 

For a wall heat rate given as qo = qo(S): The wall- 
temperature distribution is calculated by the method 


neglecting dissipation as 7) = 7(.S), and the actual 
wall-temperature distribution is 
- wise —— U? 
Vacr = TAS) +71 +37 - 
2¢/ cp 
Tre 7 U;? 
To(S) + 7;, —(1 — 7) (39) 
22 cp 


EXAMPLE 


Example results of this calculation method are shown 
for a three-dimensional source shape or half-body. The 
free-stream velocity distribution for this half-body is 
depicted on Fig. 3 with the predicted boundary-layer 
momentum and displacement thicknesses. The pre- 
dicted skin-friction and heat-transfer coefficients are 
shown on Fig. 4 in the form of the usual dimensionless 
parameters. 

As a matter of interest, the results obtained here are 
compared with similar results presented by Scholke- 
meier® for the same half-body. The comparison is 
made on the basis of the displacement thickness and 
the skin-friction coefficient. It can be observed that 
the two results are in good accord near the stagnation 
point and immediately downstream therefrom; how- 
ever, at a greater distance from the stagnation point, 
the results diverge somewhat. This divergence is the 
result of the selected velocity profiles in the boundary 
layer, Scholkemeier making use of the Pohlhausen P-4 
profiles, while the present method considers the Hartree 
“wedge flow’ profiles. Although insufficient experi- 
mental data exist to enable a selection of one velocity 
profile over the other, the use of the Hartree profiles is 
preferred by reason of the increased versatility in the 
solutions of the energy equation to account for a variety 
of surface temperature and surface heat-rate variations. 
Enough experimental evidence exists!” to substantiate 
the results of such calculations. 


A prediction of the heat-transfer coefficient for the 
half-body with a constant surface temperature is shown 
on Fig. 4 in comparison with similar results for a flat 
plate. It can be observed that the local values of the 
heat-transfer coefficients on the half-body approach 
the value for the flat plate in the limit wherein the ve- 
locity along the half-body becomes constant. 
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SUMMARY 


A calculation method has been presented which wil] 
enable a rapid prediction of the distribution of the locaj 
skin-friction and heat-transfer coefficients to be made 
for rotationally symmetrical bodies having arbitrarily 
specified free-stream velocities and surface-temperature 
variations for constant-property, incompressible, lami. 
nar boundary-layer flows.t 

Example results utilizing the proposed method haye 
been presented in the form of the boundary-layer mo- 
mentum and displacement thicknesses, the skin-friction 
coefficients, and the heat-transfer coefficients for 4 
three-dimensional half-body with constant surface 
temperature. These results utilizing the Hartree ve- 
locity profiles have been compared with similar results 
obtained by Scholkemeier utilizing the Pohlhausen P-4 
velocity profiles. 
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Statistical Aspects of Dynamic Loads 


Y. C. FUNGt 
California Institute of Technology 


ABSTRACT 


In many dynamic stress problems, the knowledge of the phys- 
ical phenomena is not precise enough to justify exact predictions 
with respect to the results of each individual observation. On 
the one hand, the loading is never exactly known, and on the 
other, the physical system may be so complicated that detailed cal- 
culations for the prediction of dynamic responses are difficult 
and the numerical results are often uncertain. In such cases the 
physical data can be regarded as a set of statistical data of ran- 
dom processes. Logical statements of the forcing functions are 
therefore statistical. The chief object of an analysis is then to 
draw valid inferences about the statistical statements of the re 
sponses. 

In this paper the statistical theory is illustrated by the gust- 
A method of determining the mean 
Through a study 


and landing-load problems. 
intensity of the dynamic response is presented. 
of the distribution of the extreme values, an 
senting the most probable largest stress (or other response quan- 
tities of interest) at any point in the structure can be found 
Such an envelope, if used as a basis for specifying the design 
stresses, can be regarded as one that yields, in the long run, a 
uniform factor of safety over the entire structure with respect 
The probability of exceeding a 


“envelope” repre- 


toa random forcing function. 
given maximum stress being known, the life expectancy of an 
aircraft can be determined as a function of the specified load in- 


tensity and vice versa 


SYMBOLS 


b = semichord of airfoil 
¢ = wing chord 
U = flight speed of airplane 


w = circular frequency, rad. per sec. 

gust speed, component of turbulent velocity normal 
to the flight path 

vertical translation of airplane 


S = wing area (sq.ft.) 

L = scale of turbulence 

k = wh/U = reduced frequency based on wing semichord 

§ = wl/U = reduced frequency based on scale of turbulence 
s = b/L, ratio of wing semichord to the scale of turbulence 
| = wingspan 


(1) INTRODUCTION 


peerage AIRCRAFT DESIGN is an art well estab- 
lished, it is still useful to re-examine periodically 
the basic point of view from which the design criteria 
This paper is concerned with one 
the dynamic loads as they 


are formulated. 
phase of aircraft design 
arise from flying through gusty air, landing impacts, 
take-off over rough surfaces, buffeting, etc. In par- 
ticular, the statistical aspects of these dynamic load- 


ing problems are considered. It will be shown that, 


Presented at the Structures Session, Annual Summer Meeting, 


L.A.S., Los Angeles, July 16-18, 1952. 
t Assistant Professor of Aeronautics. 


in certain occasions, it is logical to recognize the statis- 
tical nature of these problems and is perhaps wise to 
deviate somewhat from the current design philosophy. 


Statistical methods have been used extensively in 
theoretical physics, electrical engineering, and the 
theory of turbulence. Recently, their application to 
buffeting has been pointed out by Liepmann.' The 
present paper extends Liepmann’s ideas to some other 


problems. 


(Il) RANDOM ELEMENTS IN DYNAMIC STRESS PROBLEMS 


In a classic approach to the dynamic stress problem, 
the response to a specific forcing function is sought. 
Whenever it is permissible to linearize the governing 
equations, the mathematical principles for treating such 
problems are well known, although in aeroelasticity it 
usually leads to a set of complicated equations. 

For example, in the gust-load problem, a sharp-edged 
gust is often considered. Under a number of simplify- 
ing assumptions, it is possible to derive a ‘‘sharp- 
edged gust formula’ giving the normal acceleration 
experienced by an airplane in encountering such a gust. 
The simplifying assumptions are rarely fulfilled. Yet 
the formula is useful, for, as Rhode pointed out,* if the 
gust speed is based upon an “‘effective’’ value obtained by 
acceleration measurements made on airplanes over a 
long period of time and if the data are reduced accord- 
ing to the sharp-edged gust formula, the result will cor- 
rectly represent the actual load factor. The effective 
gust speed will be different, however, for each particu- 
lar type of airplane, because the effects of the simpli- 
fying assumptions are different for different airplane 
size, flexibility, flight speed, and type of construction. 
This approach has been found satisfactory until re- 
cently, when new models with high structural flexibility 
appear so frequently that the measurements of the 
effective gust speeds cannot keep pace with the designs. 


A theoretically more satisfactory procedure is to 
measure directly the atmospheric disturbances (gusts) 
and then design the structures accordingly. It is gen- 
erally thought that a few typical types of gust profile 
can be defined which will produce the critical dynamic 
stresses in the structures. The intensity of these typi- 
cal gusts can then be determined statistically by 
flight measurements. Thus the gust-loading problem 
is reduced to a well-defined problem of finding the 

t See reference 1 for a brief historical review and references to 


fundamental works. See also references 2-5. 
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response of an elastic airplane to some specified gust 
profiles.’ 

In applying such an analysis to the practicai design 
of an aircraft, one must know how to specify the gust 
profile. A deeper consideration will lead to the ques- 
tion whether such a formulation is sensible in physical 
reality—that is, whether such a profile can ever be 
determined by experimental measurements. 

Fig. 1 shows a record of gusts measured by an ane- 
mometer 146 m. above the ground (suspended in space 
by a balloon that is tied to the ground by three wires).f 
It is similar in nature to the anemometer records of 
wind-tunnel turbulences. From such a figure it is 
indeed impossible to say what kind of an isolated gust 
profile can characterize the real picture. Although 
1 This figure is obtained by MacCready.’ The upper curve is 
the vertical component of the gusts; the lower curve is the hori- 
zontal component. 


1953 


there is a certain regularity in each record, it is impos- 
sible to predict the exact shape, height, and duration of 
-ach peak. 

A similar question arises in other dynamic stress 
problems. In practice, most dynamic stress problems 
are characterized by the unpredictable nature of the 
forcing functions. The exact form of the forcing fune- 
tion is always uncertain to some extent because of the 
large number of variables entering the problem. Thus, 
when a few important ones are considered in a theory, 
other variables appear as external disturbing influ. 
ences. The random nature of the atmospheric turbu- 
lence as shown in the records of Fig. | is due to the in- 
fluences of the viscosity, pressure, density, temper- 
ature, humidity, and velocity distributions, varying in 
such a complicated manner that mechanistic predictions 
based on the hydrodynamic equations are practically 
impossible. 

Consider the landing impact as another example. 
In addition to the landing speed, drop height, and the 
airplane attitude as ordinarily considered in the landing 
analysis, there are a number of factors beyond the con- 
trol of an airplane designer, such as the softness of the 
ground, the tire pressure, the change of temperature 
and its effect on the viscosity of the fluid in the oleo 
strut and the air cushion, the surface friction, the 
bumpiness of the ground, the pilot technique, etc. 
The list of possible factors is almost unlimited, and it 
is of little interest to compile such a list because it 
would not be useful. In the landing analysis, the as- 
suinption usually made is that it is possible to introduce 
average values of these factors and that the practical 
variations of these factors from their average values 
introduce only small errors. When 
stresses due to the variations of the neglected variables 


the secondary 


are not negligible, we are left with two ways to account 
for them: (1) by introducing a sufficiently large factor 
of safety (or a higher load factor), or (2) by studying the 
variations and introducing some allowances based on 
statistical inferences. The second alternative is, of 
course, more attractive if it can be accomplished with- 
out causing too much difficulty in the analysis. 

In Fig. 2 are shown the experimental results of Ram- 
berg and McPherson of the National Bureau of Stand- 
In their tests a 
model single-spar wing having only flexural degrees of 
The model was 
The land- 
ing conditions were varied by changing the softness of 
the rubber on which the model was dropped, the initial 


ards on the landing impact force.® 


freedom was fitted with an oleo strut. 
dropped from a carefully controlled height. 


spring loading, and the opening of pressure relief ports 
in the damper. The transient bending stresses were 
measured both at the wing root and at an outboard 
station. The test results are summarized in Table | 
and Fig. 2. 


in the maximum impact force, bending moment, ac- 


It is seen that considerable variations exist 


celeration at the fuselage, and the shape of the impact 


force curves. These variations will cause variations 
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TABLE | 


Maximum Values of the Applied Impact Force (Sum of Spring and Damper Forces), Bending Moments at 1.5- and 14.5-In. Stations, 
and Accelerations at Root (froii Reference 9) 


Maximum 

Impact Bending 

Force per Moment 
Half-Wing Mi.1.5 

Record (Lbs. ) (Lbs.-In.) 
a 95 750 
b 105 1,130 
c 100 1,060 
d 92 970 
€ 120 1,370 
f 113 1,170 


in the dynamic stresses in the airplane, especially if 
the airplane structure is flexible. Hence, if a typical 
impact force-time curve is chosen among all the curves 
of Fig. 2, it is necessary to increase the amplitude of 
this curve by a factor (of safety) of considerable magni- 
tude in order to overshadow the maximum stresses that 
may be induced by other impact force-time curves. 
Under such circumstances it seems desirable to make a 
more rational approach so that a more realistic design 
specification may be obtained. It is from consider- 
ations of this nature that a statistical theory is origi- 
nated. 


Let the curve shown on the left-hand side of Fig. 3 
be taken as a typical average curve. The individual 
impact curves of Fig. 2 can be represented by super- 
positions as shown in Fig. 3. An essential feature of 
the curves on the right-hand side of this figure is their 
irregularity. Remembering that irregularities 
are caused by a large number of factors, it is not sur- 
prising that no definite rule can be stated which will 


such 


predict their possible form. 


If the uncertainty of predicting the forcing function 
is recognized, the problem becomes statistical. One 
then tries to state the main features of the response, as 
well as those of the forcing function, in terms of statisti- 
cal averages and probability distributions. 


(111) THE PROBABILITY DISTRIBUTION 


It was stated in the last section that in some dynamic 
stress problems our knowledge of the physical phe- 
nomena is not precise enough to justify exact predica- 
tions with respect to the results of each individual ob- 
servation. On the one hand, the loading is never ex- 
actly known, and on the other, the physical system is 
so complicated that detail calculations for the predic- 
tion of dynamic responses are difficult and the numeri- 
cal results are often uncertain. In such cases, we shall 
regard the physical data as a set of statistical data of 
random variables. The chief object of an analysis is 
then to draw valid inferences from the data. 

t This curve is the mean of all the curves of Fig. 2—i.e., the 
mean value of the impact force at each instant of time following 


the contact of the landing gear with the ground. 


Bending Maximum 
Moment Acceleration Impact 
Mo, 14.5 at Root Duration 
(Lbs.-In.) (In./Sec.?) g (Sec. ) 
370 1,580 4.1 0.086 
650 1,930 5.0 0.076 
550 1,580 4.1 0.062 
540 1,510 3.9 0.060 
820 2,220 5.8 0.052 
600 1 ,620 4.2 0.050 


The statistical nature of a random variable is char- 
acterized by its distribution function.. For a random 
variable y, the distribution function F(x) is the prob- 
ability of the event y < x: 

F(x) = P(y < x) (1) 


Obviously, F(x) is a nondecreasing function, and 


F(—©) =0, F(+ ”) =1 (2) 


If the derivative F’(x) = f(x) exists, as will be assumed 
hereafter, f(x) is called the probability density or the 
frequency function of the random variable. A prob- 
ability density is nonnegative and has the integral 
l over (— ©, ©), since, by definition, 


Pia<y <b) = F(b) — F(a) (3) 


In the limit it is seen that the probability that the 
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NORMAL PROBABILITY DENSITY f(x) 
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Normal distribution. 


Fic. 4. 


variable y takes a value belonging to the interval 
xs<y<ne + Ax 


is for small Av asymptotically equal to f(x) Ax, which is 
written in the usual differential notation 


Pix <y<x + dx) = f(x) dx (4) 


Moments 


The integrals 


a, = f a ie) de; v= 1, 2).8,.5...) (5) 


if absolutely convergent, are called the first, second, 
third, ..., moment of the distribution function. The 
first moment, called the mean, is often denoted by the 


letter m: 
m = / xf(x) dx (6) 
The integrals 


by = / (x — m)” f(x) dx (7) 


are called the central moments. 

The mean m is a kind of measure of the location of the 
variable y. It gives a typical value of the variable, 
which may be conceived as a central point of the dis- 
tribution. 

The second central moment pe gives an idea of how 
widely the values of the variable are spread on either 
This is called the variance of the 
always pw > O; 


side of the mean. 
variable. We have, of 
ue = 0 if, and only if, the variable assumes only one 
value. Generally, the smaller the variance, the greater 
the concentration, and conversely. 

In order to obtain a quantity of the first dimension 
in units of the variable, it is often preferable to use the 
nonnegative square root of w2, which is called the stand- 


course, 


ard deviation oi the variable and is denoted by a, 


9 


o? = wo = ag — m? (8) 
Example: The Normal Distribution 


The normal distribution function is defined by the 


relation 


1953 


1 x 
Wa) — os ~13/2 
F(x) = ss I. e dt (9) 


The corresponding normal probability density is 


f(x) = (1/V/2n) e7*"”” (10) 


Diagrams of the function f(x) is given in Fig. 4. The 
mean value of the distribution is 0, and the standard 
deviation is 1. 

A random variable & is said to be normally dis- 
tributed with the parameters m and o—or in brief, nor- 
mal (m, o)—if the distribution function of & is F[(x - 
m)/o|, where F is given by Eq. (11) and o > 0 and m 
are the standard deviation and the mean, respectively, 
The probability density is then 


f(x) a (1 ov/ 2r) e © m)?/2a? (11) 


(IV) STATISTICAL AVERAGEST 


Thus, 
the mean and the mean square of a function of time 


The concept of time averages is well known. 


u(t) over a time interval 27 are 


1 *to+7 
Gi = - ff , ¥O dt (12) 
| “to +1 
(u — u)? = =f , u(t) dt — nu” (13) 
— h- 


If a and (w — %)* are dependent neither on f) nor on T, 
provided that 7 is sufficiently large, then u(t) is called 
a stationary random process. 

A set 
In the 
example of landing impact measurements such as 


The concept of ensemble average is important. 
of observations forms an ensemble of events. 


those represented in Fig. 2, each impact force-time 
curve is a member in an ensemble of such curves. 

Imagine that a large number of experiments be made 
simultaneously on identical models under similar con- 
ditions. Let the records be numbered and denoted by 
U,(t), U2(t),..., Uv(t). If the total number of records is 
N, the following averages may be formed: 


t= (uw + w+... + uy)/N (14) 
N 
(u—#)?= D> (u; — a)?/N (15) 
i=1 
Assume that # and (u — “#)? tend to definite limits m 


and o’, respectively, as V— ©. Then a, (u — m)’, or 
m, o* are ensemble averages of the random process u(/). 

If the ensemble averages of u(t) are independent of 
time, then u(t) is said to be stationary. For a stationary 
random process, the ensemble averages are identical 
with the time averages. 

There are many flows of interest that are approxi- 
mately homogeneous and stationary when the time 


t See references 1-5 for greater detail and more precise dis- 


cussions. 
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and distance scales are properly chosen. Such is the 
case of atmospheric turbulence within a suitable ex- 
panse of time and space. 

It is natural to expect that the behavior of a dynamic 
system depends on the manner in which the excitation 
changes in successive instants of time. Simply know- 
ing the mean and intensity of the excitation is not 
sufficient to establish the corresponding mean and in- 
tensity of the response. Additional information ex- 
pressing the sequential time characteristics of the ex- 
citation is necessary. One of the measures describing 
the sequential behavior has its origin in the harmonic 
analysis and is called the power spectrum. Alterna- 
tively, a more easily defined quantity called the cor- 
relation function may be used. 


The Power Spectrum 

Let us consider a stationary random process u(/) with 
mean m = O and define the mean square of f(/) over 
an interval 27 as in Eq. (13). 

Let u(t) be truncated so that it becomes zero outside 
the interval (—7, 7). Let the truncated function be 
written as u7(t). The Fourier integral of u7(¢) exists, 


and we may write 


a 

ur(t) = 7] | F(w)e dw (16) 
Vv or x 

F we : tool (17) 

F(w) \/ on | , ler (t) e dt é 


It is well known (Parseval theorem) that 
yx 7 on: 

/ ur*(t) dt = / ur (t) dt = / F(w)|? dw 
in 7 en 


Hence, 


f- F (w) |?d fr He) de (19 
u(t) = =. Js (w) |\-dw = a= dw ()) 
27 ~ 0 7 


For a stationary random process u“*(f) tends to a con- 
stant, so |F(w)|*/7 tends to a definite limit. Let 


(18) 


F(w) |? 
p(w) = lim = (20) 
sag 7 


u*(t) = y p(w) dw (21) 
0 


Following Taylor, p(w) is called the power spectrum of 
u(t). The element p(w)dw gives the contribution to 


then 


u*(t) from simple harmonic elements with frequencies 
ranging from w to w + dw. The integral 


| P(w) dw 
0 


represents the contribution to u°(f) from frequencies 


less than w. 


Correlation Functions 


Consider again a stationary random process u(f). 
The average value 


| 7 
V(r) = u(t)u(t + r) = af u(t) u(t + 7) dt (22) 


=f 
is well defined, provided that 7 is sufficiently large. 
¥(r) is a function of 7 and is called a correlation function. 


The following properties of the correlation function 
can be easily seen: 


¥(O) = u(t) 23) 
lim ¥(7) = 0 (if a = O) (24) 
¥(r) = ¥(—7) (25) 


Moreover, the reciprocal relations between the power 
spectrum and the correlation function can be estab- 


lished! 


V(r) = I p(w) cos wr dw (26) 
J0 
; f° ‘s 
p(w) = ¥(r) cos wr dw (27) 
wr J0 


(V) THe POWER SPECTRA OF THE EXCITATION AND 
THE RESPONSE 


Suppose that the response y and the excitation f are 
connected by the differential equation 
If.f(t) is a harmonic function 
f(t) = Ae (29) 
where A is real, the solution y(‘) which represents a 
steady motion is 
y(t) = Ae /Z(iw) (30) 
where 
Z(iw) = wo”? — w*® + Bw (31) 


is called the impedance of the system represented by Eq. 
(28). Generalizing this result to an arbitrary forcing 
function f(t) represented by a Fourier integral 


twl *2*) 
ff = V3 A(w) e dw (32) 
<7 x 


we have the response 


“2 A(w 
. =: e dw (33) 
V 24 o Z(tw) 


which is valid if the integral converges. 


y(t) = 


If f(t) is a stationary random function, it does not 
tend to zero as t > © and its Fourier Transform does 


not exist. But a truncated function can be defined so 
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that Eqs. (32) and (33) hold within a large time in- 
terval (—7, 7). The following relation can be proved: 


~ Dp 
y*(t) = f l f (w) a dw (34) 
: 0 |Z(tw)|? 


The phase 





where p(w) is the power spectrum of f(f). 
relationship between the excitation and the response, 
represented by the argument of the complex number 
Z(tw), is completely obliterated in the power spectra 
relations. 

The preceding conclusion holds for a much wider class 
of dynamic systems than that represented by Eq. (28). 
It holds for higher order linear differential equations, 
linear integral equations, or linear integral-differential 
equations under mild restrictions. 


(VI) THE RESPONSE OF A RIGID AIRPLANE TO 
ATMOSPHERIC TURBULENCES 


As an illustration of the statistical approach, let us 
consider the motion of an airplane in response to atmos- 
pheric turbulences. For simplicity, let us introduce 
the assumptions that the airplane can be regarded as 
rigid, having uniform forward speed U’ and responding 
only to gusts normal to the flight path, and that only 
the vertical translation degree of freedom is permitted. 

The atmospheric turbulence (gusts) may depend on 
the condition of flow and the thermal stability of the 
air. Therefore the statement of the characteristics of 
turbulences must be based on measurements relevant 
to specific meteorological conditions.t However, Clem- 
entson'® has shown that the correlation functions 
(and the power specra) of atmospheric turbulence in 
different conditions (unstable air mass, water-land 
discontinuities, thunderstorms, and mountainous ter- 
rain) have remarkable similarity to each other. 

The following analysis will not refer to any particular 
weather condition but to the isotropic turbulence as 
measured in wind tunnels. Obviously, however, the 
analysis can be extended to other turbulence fields. 

It is assumed, therefore, that the gust field has a 
power spectrum given by! 


L 1+ 32 


wl 
p(w) = w? oe 


2 (35) 
awl (1 + £7)? U 


where w’ is the mean square (intensity squared) of the 
gust, L is the scale of turbulence, and UL’ is the speed of 
the general flow. 

The scale of turbulence is defined by the area under 
the correlation curve, 


L= f° Roar (36) 
0 


+ At the present time no general conclusion has been reached in 
relating the distribution function of the gust velocity components, 
the power spectra, and the correlation functions to the velocity 
gradient, lapse rate, and other meteorological factors. For avail- 
able information, see, for example, references 10-12 


where 


w(x, t) w(x + 7, t) 


R(r) = (37) 


w 

Let us consider a strip of unit length of an airfoil of 
infinite span in a two-dimensional flow. The mean- 
flow velocity is U parallel to the mean position of the 
airfoil chord. The thickness and the amplitude of the 
displacement of the airfoil are assumed small in com- 
parison with the chord. To the mean velocity of flow 
U is superposed a small turbulent motion that has only 
one component of velocity w in a direction normal to 
the chord. The turbulent motion in the fluid is trans- 
ported with the mean flow, so that w is given by a func- 


tion 
w= wi{t — (x/U)] 


The lift of the airfoil is partly induced by the disturb- 
ance w and partly by the motion of the airfoil itself. 
Within the framework of the linearized theory, the lift 
induced by these two parts are superposable. The 
equation of motion of the airfoil (having the transla- 
tional degree of freedom only) can be written as 


mz + Liz) = —L(w) (38) 


where ,L(w) indicates the lift corresponding to w(t) 
and L(z) that for the motion of the airfoil. It is clear 
that the problem of the response of an airfoil to a turbu- 
lent flow can be separated into two problems: (1) the 
determination of the lift produced by the turbulent 
flow on an airfoil which is stationary relative to the 
mean-flow velocity, and (2) the determination of the 
disturbed motion of the airfoil to an exciting force, with 
the flow regarded as uniform and without turbulence. 

Let p,(w) and p(w) denote the power spectra of the 
turbulent flow and the lift, and let x,(w), x;(w) denote 
the square of the absolute values of the frequency re- 
sponse (admittance) of the lift and the airplane acceler- 
ation, respectively; then 


Pilw) = xalw) pp(w) (39) 
Pw) = xs(w)xa(w) py(w) (40 
Hence the problem is really the determination of the 
quantities x,(w) and x,(w). 
Let the vertical gust velocity be 
w(x, t) = wy b@!-2/) (41) 
The expression of lift is given by Sears! 
L = 2rpbU w e™ o(k) (42 
where k is the reduced frequency 
k = wb/U (43) 
and ¢(k) is a complex function 


e(k) = [Jo(k) — 1Ji(k)] C(k) +i Si(k) (44 


Ji, J; are the Bessel functions of the first kind and 
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C(k) is the Theodorsen function. The resultant lift 
acts through the one-quarter-chord point from the 


leading edge. The factor 2rpbl’ g(k) represents the 


frequency response (admittance) of the lift. 
Liepmann has shown! that a simple approximation 


of | ¢(R)|? is 


o(k)|? = 1/(1 + 2nk) (45) 


Admittance of the Airplane 
The equation of motion of a two-dimensional airfoil 
subject to a sinusoidal force can be reduced into the 


following form: 


U? . apU? S : : ~ 
. b? k? + 9 bh [2C(R) ik _ k*)¢ Zo = Po (46) 
where P = Pye’ is the exciting force and z = ze" is 


The quantity in the braces { } 


Its inverse 


the wing displacement. 
is the complex impedance from F(t) to 2(f). 
multiplied by —w? = —(Uk/b)? gives the admittance 
of the acceleration z. Let 


k = 2m/mpSb (47) 


be the density ratio of the airplane; then the admittance 
can be written as 


Z 2 k 


—— : (48) 
P wpSb (1 + «wk + 2G — 2F 


The Gust Response 


The intensity of the acceleration V 3° can now be 
Assume the power spectrum of the gust 
to be given by Eq. (35). Neglecting the lift on the 
tail surfaces, and using strip theory, the lift on the wing 
due to a gust can be obtained by multiplying Eq. (42) 


computed. 
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by the wing span. Then in the notation of Eq. (40), we 


have 


Xe = (wpSU)?| y(k)|? (49) 


5 |2 9 \2 
,=f/=| =(— x 
os r (5) 


k? 


50 
4(F? + G2?) + 4GRK1 +e +R +42 


The mean square of the acceleration is, therefore, 


- wh 
2 = / Xa(R)xXs(R) py (w)dw, k= v (51) 
0 


The exact integration of the above expression is diffi- 
cult because of the complicated manner in which the 
Bessel functions are involved. To obtain an approxti- 
mate solution, the approximation for |¢(k)|*? given by 
As to xs, we note that, in the 
Hence 
More- 


Eq. (45) may be used. 
full range (0, ©), G? is much less than F°. 
we may neglect terms involving G? in Eq. (16). 
over, since 


|F| < 1, G| < 0.2 


' 
while in practice, x varies from a number of order 40 
for trainers and transports to 150 for high-speed fight- 
ers, a fair approximation of Eq. (16) is simply 


2 . k? 
X: = ( ) (52 
apSb t+ k7(1 + x)? 


Using the approximations, Eqs. (45) and (52) and intro- 
ducing the notations 


s = b/L (53) 
a = 2/s(1 + x) (54) 


where s is the ratio of the wing semichord to the turbu- 


lence scale, we obtain 


4 U/? “~ g l 1 + 3é _ 
32 = yw? = é - = dé (O09) 
wb*(1 + x)? Jo a> + & 1 + 2asté (1 + €)? 


By partial fractions, 


ge 1+ 3¢ 


+ 
e+eate? a+82 142° (+2)? 


where 


a*(3a* — 1) 
A= = —) B 
(a* — 1)* 


Furthermore, 


] « dé 
(a, s) = ~~ = — 
; o (1 + 2ast) (a? + £) 


ja? — 3 2 (56 
: ae C=-7, 00) 
(a* — 1)* a? — | 
2ars l ) > 
» 2 y (O¢ 
1 + 4n°s*a’ log 2es + loge + las 


Differentiating both sides with respect to a and then dividing through by (— 2a), we obtain 








Srs3 


yf dé 
Ip(a, s) = de he ae 
o (1 + 2xsk) (a? + €)? (1 + 47r?s?a’) 


Hence, from Eq. (55), 


32 = w? (4U"/rb*(1 + «)?] J (a, s) (59) 


where J(a, s) = Alj(a, s) + B(1, s) + Cha(1,s) (60) 
,(1, s), J2(1, s) are given, respectively, by Eqs. (57) and 
(58), with the parameter a replaced by 1; and A, B, and 
C are given by Eq. (56). The variation of /(a, s) is 
Note that 2? ~ 0 both when s > 0 

Hence, as the scale of turbulence 


shown in Fig. 5. 
and when s > ©. 
becomes either negligibly small or infinitely large in 
comparison with the wing chord, the intensity of the 
acceleration experienced by the airplane will tend to 
zero. This can be expected, because when the wing 
chord is extremely large in comparison with the turbu- 
lence scale, the “‘gusts’’ are smoothed out by canceling 
each other over the wing. On the other hand, when 
the chord length is extremely small in comparison with 
the scale of turbulence, the airfoil behaves quasi- 
stationary. Since a rigid airplane having only the 
translation degree of freedom will experience no ac- 
celeration in a stationary flow, the limiting case 22 > 
0 as s > Ois to be expected. 


l 
9 . y “ 
(tog 2rs + loga + x) + 


TS ( 1 
(1 + 42*s*a?)a? \4sa 


o 1) (58) 


The factor J(a, s) shows clearly the effect of the size 
of the airplane relative to the scale of turbulence. The 
effect of finite span is neglected. Liepmann! has shown 
that essentially the finite span effect is to multiply 
the right-hand side of Eq. (59) by a factor 


(L/l) [1 — e7”*] 


(61) 


where / is the length of the wing span. 


(VII) HOuUSNER’S “STRESS ENVELOPE”’ 


In structural designs, it is often required to find the 
maximum stresses expected at any point in a structure, 
In this connection it is interesting to notice an approach 
used by Housner™ in the study of strong-motion earth- 
quakes, which is a nonstationary random process. 

Consider the following specific problem: A beam is 
clamped at one end and free at the other. The clamp- 
ing wall is subjected to a variable acceleration that is a 
random function of time. Find the response of the 
beam. 

If the variable acceleration of the wall is a(t) and if 
a(t) = Ofor?t < O, it is not difficult to derive the beam 
displacement u(x, f), 
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K,; ' . 
u(x, t) = YB f(x) a(r) sin w(t — 7) dr (62) 
‘ 5 wy 0 
where 
f(x) = the 7th normal mode of the beam 
wy = the circular frequency of vibration of the 7th mode 
K, = JS fidp/S f2 dp, dp being a mass element of the beam, and the integrals are taken over the entire mass 


This can be written as 


K, nt t 
u= > —fi(x) | sin wit y a(r) cos wit dt — COS wit a(r) sin w;t dr 
Wj 0 0 


1 
(63) 
K, j “t 2 t 2 : 
= > —f,(x) ) : a(r) cos wit dr} + a(r) sin w;7 dr ri sin w; (t — aj) 
i 0 0 
where a; is a phase factor that depends on the values of the integrals. Introducing the notation 
t 2 t 27)'2 
Rg, = ( a(r) COS w;T dr) + ( a(r) sin wT dr) (64) 
0 0 
Eq. (63) becomes 
kK; - 
ilies > fi(x) Rj sin w(t — aj) (60) 
i Wi 
The quantity R; is a function of /, a(/), and w;. It displacement at any point is then 
characterizes the acceleration a(/). K 
Housner'® ' has shown that, if a(¢) represents the u=R>Y—|fix)| (67) 
i 


ground acceleration encountered in strong-motion 
earthquakes, the maximum values of R; (reached at 
certain time /) vary erratically with the frequency «;. 
For two frequencies close to each other, R; may have 
markedly different values. On the basis of analyses of 
many past earthquake records, it cannot be concluded 
that the maximum R; will necessarily be larger at one 
frequency of vibration than at another. In fact, on 
averaging a number of such R; vs. w; curves, it is found 
that the maximum R; is essentially independent of w;, 
except for very high frequencies (R; ~ 0 as w; > ©). 
Housner shows" that this should really be expected if 
the earthquake waves consist of a series of impulses 
arriving at random intervals. 

If R; is a constant independent of w;, we may write 
R, as R, which is a measure of the intensity of the earth- 
quake. In this case Eq. (65) may be rewritten in terms 
of R, 


K, a 
u=R>D f(x) sin w(t — a;) (66) 
1 Wj 


Examining the trigonometric term, it is seen that, be 
cause of the variable phase factor a;, the effect of the 
modes may be additive or they may tend to cancel each 
other. Either of these possibilities may occur, and 
whether or not they do depends upon the earthquake. 
The worst case is when the maximum contributions of 
that is, when sin w;(¢ — a;) = | 
The maximum 


each mode are additive 
and f(x) has its absolute value | f;(x)). 


It should be noted that Eq. (67) does not represent the 
actual displacement at any time, but at any ¢ it repre- 
sents the maximum that would occur if all modes 
were additive. In this sense, Eq. (67) may be said to 
represent the envelope of maximum displacements. 

The argument that a certain envelope can be con- 
structed which represents the worst possible combina- 
tion of all the modes can be applied to the more general 
case in which R; does depend on the frequency «,. 
Then the envelope of displacement at any point is 
K 

| Ri(w,)| | fi(x) (68) 


W; 


u= > 


F 


Example: A Uniform*Cantilever Beam 


Consider a uniform cantilever beam that has a uni 
form distribution of mass and rigidity and which de 
flects with bending deformation only. Let the mass 
per unit length be m and the flexural rigidity be //. 
The free vibration modes f;(x) and the frequencies w, 
are well known and will not be quoted here. 

The envelope of maximum displacements for a canti 
lever beam is given by Eq. (67) if we assume random 
impulses so that R is a constant. The envelope of 
maximum bending moments is given by 

d*u K; | d*f; 
M = El = REID “7 RI V mEI Xy, 
d 5; Ww dx 


(69) 
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where 
K,| df, , 
x m a fa (70) 
rage Midd Ot 


and /p; is the 7th root of the equation 
cos p:cosh p + 1 = 0 (71) 


The shape of the envelope of maximum bending mo- 
ments is given by X,,. The absolute magnitude de- 
pends upon the intensity of disturbance R and upon the 
mass per unit length and stiffness, V mE/. The shape 
X, is shown in Fig. 6 for the first mode, the sum of 
the first three modes, and the sum of the first six 
modes. 

On the same figure, for comparison, the shape of a 
bending moment curve corresponding to a uniformly 
distributed load is shown by the dotted line. 


Similarly, the envelope of maximum shearing force 
can be obtained from the equation 
d*u 


_ _ K;, | df, 
ger SD, las 


(72) 


S = El 


The shape of the maximum shear envelope is given in 
Fig. 7. The shear that would have developed if a uni- 
formly distributed load is applied statically is shown 
by the dotted line in the same figure. 


(VIII) ExTREME VALUEST 


In a sample of » values of any random variable, there 
are always two finite and uniquely determined extreme 
values. Let us consider the case when the population 
has a distribution of the continuous type, with a dis 
tribution function F(x) and a_ probability density 
f(x) = F’(x). Let x denote the largest value in a sample 
of nm from this population. 
g(x)dx in the sampling distribution of x is identical with 
the probability that, among the » sample values, 
n — lare < xand | falls between x and x + dx. 
by the rules of multiplication of independent combined 
events and the definition of the distribution function, 


The probability element 


Hence, 


we have 


g(x) dx = n [F(x)]"~' f(x) dx (73) 
Let us introduce a new variable by the substitution 
&=n [1 — F(z)] (74) 
Since F(.) varies between 0 and 1, we have 
O<ét<n (75) 
By differentiation, di = —n f(x) dx. The expression 
(73) is therefore transformed into 
g(x) dx = —[1 — (&/n)]"~' dé (76) 


- 


The probability density of the new variable é is, there 


fore, 

h(é) = [1 — (&/n)]"~' (77) 
for 0 < & < n and h(E) = O outside (0, m). Asn 
o, h(£) converges for any  < 0 to the limit 

lim h(é) = e~* (78) 


n-—>« 


When the initial distribution function F is given, it is 
sometimes possible to solve Eq. (74) explicitly for x. 
The extreme value x is then obtained in terms of an 
auxiliary variable £ whose probability density is known. 
Generally, however, such explicit solution cannot be 
given, but the limiting form of the distribution of the 
largest value when 7 is extremely large can often be 
found. 

For example, if the population distribution is normal 
(m, o) and if x is the largest value in a sample of n values 
from this distribution, Eq. (74) gives 


n i , = 
— —(t — m)2/2e2 ( 
; = , oe ee (79) 
oV 2r x 


The asymptotic solution of this equation, for large 


values of n, is 


+ The method that follows here is given by Cramer, reference 
18, pp. 370-378, where further references to papers by Fisher, 
Tippet, Irwin, Pearson, Davies, and Hartley are given. Furthet 
references to a series of papers by Gumbel, Kimball, and Kendall 


can be found in references 19 and 20 





wh 


Th 
the 


wh 


Th 


is t 


def 


for 


me 


an) 
Hei 


whe 
fun 


App 
I 


vel 
cele 
prol 
cial 
terr 
to ¢ 
spec 

F 
will 


Idet 


As s 
and 


here 


reme 
tion 
dis- 
sity 
uple 
lent 
vith 
lies, 
nce, 
ned 
ion, 











STATISTICAL ASPECTS OF DYNAMIC LOADS 


log log log 4 
; og og n + og 74 


2 V2 log n 
Co 


Vv 
V2 log n 


y=mta V2 log n — 


(SO) 


where 


v = —logé (S1) 


, = 


The variable v, by virtue of Eqs. (Sl) and (78), has 
the probability density 7(v) : 


jv) =e" * (82) 


which is shown in Fig. 8. The distribution function is 


J(v) =e ° (83) 


7 = 


The complement /*(v), which gives the probability that 


v be larger than a specified value, 


J*(v) =1—e° | (84) 


is tabulated by Press in Table | of reference 19. 
As another example, consider the Laplace distribution 
defined by the probability density 


f(x) = (1/2A)e— *-*"" (S85) 


for which the mean is yp, the first moment about the 
mean is zero, and the second moment about the mean 
is 24°. The distribution function F(x) is obtained by 
an integration. Eq. (74) gives, when x > u, 
E=n {1 — F(x))] = @/2e *™ . (S6) 


Hence, 
x = wp +A log (m/2) + Av (S7) 


where v = —log é has exactly the same distribution 


function as given by Eq. (83). 


Application to Gust-Load Problem 
It was shown by Press’ that the “‘effective’’ gust 
velocity and the peaks of the increment of normal ac- 
celeration experienced by an airplane in flight have a 
probability distribution that is approximately Lapla- 
cian. Therefore, if the parameters u and ) were de- 
termined from flight-test records, it would be possible 
to calculate the probability of encountering a gust of 
specified intensity. 
Following Gumbel and Press, the transformation (87) 
will be written as 
v= a(x — u) (88) 
Identifying with Eq. (87), it is seen that 
a=1/dr (S9) 
u = p+ d log (n/2) (90) 


As shown in Fig. 8, the frequency function j(v) has one, 
and only one, maximum, reached when v = 0. Hence, 


tw 
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Fic. 8 The probability density 7(v) 


v = Ois the so-called mode of the distribution 7(v) and 
is in a sense a measure of the most probable value of 
v. From Eq. (S&) it is seen that v = 0 implies x = u, 
and, hence, « represents a most probable value of the 
extreme value x in a sample of size n. 

In practical applications to the gust problem, samples 
of the extreme values are readily obtainable from the 
N.A.C.A. V-G Recorders. From a V-G record it 1s 
impossible to tell the time history of the gusts, but the 
extreme values of the acceleration are obtained. If a 
number of such records are available, each representing 
approximately the same number of hours of flight, the 
statistical distribution of the extreme values can be 
estimated. The probability of experiencing an acceler- 
ation exceeding a specific value can then be calcu 
lated. From the so-called 
“effective” gust velocity can be calculated on the basis 


acceleration data the 


of the sharp-edged gust formula. 

In order to utilize the type of data named above, it 
is necessary to estimate the parameters a and u directly 
from the samples of extreme values. The parameters 
a and u can be determined from the mean m and the 
standard deviation ¢ of the distribution of the largest 
value from the following equations, which are given by 
Gumbel :”° 

u =m — (y/a) (91) 
l/a = (VW6/n) (92) 


where y is Euler’s constant and is approximately equal 
to 0.5772. If Eqs. (91) and (92) are assumed to be 
true when the parameter values m and @ are replaced 


by their respective sample estimates ¥ and s, 


/ 
s = VX(x — £)2/(N — 1) (93) 


where JN is the total number of observations and > 


sums over the N observed values of x, then, 
u= xX — (y/a) (94) 
l/a = (V6/n)s (95) 


For samples of data, estimates of the values of u and a 
may be determined from Eqs. (94) and (95). By 
means of the transformation law, Eq. (SS), the equiva- 
lent values between x and v are found. Since the dis 
tribution function of v is known, the probabilities that 
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a value of v and the associated value of x will be exceeded 
may then be obtained simply from the tables of J*(v). 

In an example given by Press,’ gust velocity meas- 
urements from 485 traverses of thunderstorms from 
the 1946 operations of the U.S. Weather Bureau 
Thunderstorm Project are analyzed. The parameters 
u and a are calculated according to Eqs. (94) and (95). 
The transformation (88) now appears as 


v = 0.2072 (w. — 12.837) (96) 


where w, is the “‘effective’’ gust velocity. For a given 
value of w,, the equivalent value v can be computed 
from Eq. (96), and then the associated probability of 
exceeding w, in each traverse can be found. 

From the known probability of exceeding a specified 
gust velocity, the average miles of flight required to 
exceed that gust velocity can be determined by the equa- 
tion 


mileage = (U/P)r (97) 


where L’ is the average flight speed, P is the probability 
that a given value of w, will be exceeded on each record, 
and 7 is the average flight-hours per record. 


(IX) DyNAmIc STRESSES IN AN ELASTIC STRUCTURE- 
STATIONARY RANDOM PROCESS 


For a stationary random process the time average is 
the same as the ensemble average. This remark, to- 
gether with a result on the distribution of extreme val- 
ues, gives a new significance to the root mean square of 
the random process. Eventually, it will make possible 
the determination of not only the most probable value 
of the largest stress at any point in a structure but also 
the life expectancy of a given structure with respect to 
the random process considered. 

Let us consider the gust problem. 
turbulence will be regarded as a stationary random proc- 
ess with a mean zero. The root mean square (aver- 
aged over time) of the turbulent velocity component w, 
normal to the flight path of an aircraft, is also the stand- 
ard deviation at any time of a large number of records 
ef velocity measurements performed simultaneously 


The atmospheric 


[r(x )[? = X lam)? gm2(x) + OD 


men 


the second summation being taken once over all pairs of 
(m, n) under the stipulation m # n. 

The variance of the stress component at Xis 7*. As- 
suming normal distribution, the probability density of 
the stress 7 at X is 


f(r) = (1/oV/ ale?" (103) 


where 


(104) 


o = 7(2, t)* 


1953 


under similar conditions. Similar interpretation ap- 
plies to the mean and standard deviation of the air. 
plane response (the elastic displacement, the acceler. 
ation at c.g., or the dynamic stresses; etc.). 

However, the probabilistic structure of the gust 
response is not entirely given by the mean and standard 
deviation. Additional information is required to char- 
acterize the gust response as a random process. 


It is generally assumed that the velocity component 
in a turbulent flow is normally distributed. If the 
aeroelastic system is assumed linear, the response will 
also be normally distributed. 

Let us assume, then, that the probability distribution 


The mean and stand- 
ard deviation being known, the probabilistic structure 


of the gust response is normal. 


In particular, the distribution 
of the extreme values is defined. 


is completely defined. 


Let the displacement u of the elastic body be de- 
scribed by the generalized coordinates f,(x) and gen- 
eralized displacements q,,(t), 


u (x,t) = D> gn(t) fa(X) (98) 
n 
The stress at any point in the elastic body is a linear 


function of the derivatives of u with respect to the posi- 

tion vector x. Hence, we may write 

tr (X, t) = DO gn(t) gn(X) (99) 
n 


where 7 denotes a component of stress at X at time ¢ and 


£,(X) being functions of x derived from f,(x). Accord- 
ingly, 

r(x, t)|2 = SY anlt)qn* (t) gm(X) gn(X) (100) 
where q,,*(/) is the complex conjugate of q,(¢). Taking 


averages over a sufficiently long period of time, we ob- 
tain 


i7(x, 2/2? = SO DY qmlt) Qn*(t) m(X) gn(X) (101) 


m n 


The left-hand side is a function of x. The right-hand 


side is real, for Eq. (101) may be written as 


[am(t)Qn* (t) + Qn*(t)qn(t) | Zm(X) En(X) (102) 


The maximum stress expected at any point in a 
structure can be found by the method described in 
Section (VIII). It is shown there that, if the dis- 
tribution of a random variable x is normal with mean 
m = O and standard deviation o, the maximum value 
of x in m sample values has the following probability 
density asu > o: 


(105) 


jl). 
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where 
v=a(x — u) (106) 
log log n + log 4x 
u=o V2 log n — . —— (107) 
2/2 log n 
ja = «/V2 log n (108) 
The point v = 0 is the mode of the frequency function 


i(v). Hence, the ‘“‘most probable’ value of the largest 


e 
value of x in m samples (n > ©) is, from Eq. (106), 


x=4u (109) 


Let the random variable x stand for the stress r. 
It is seen that for a fixed large value of n the most 
probable value of the largest stress is 


tmar. (%) = const. V r(x, #)? (110) 


Therefore the most probable largest stress is propor- 
tional to the root mean square of 7(x, ¢) and is a func- 
tion of the space coordinates x. 

A surface representing the most probable largest 
stress as a function of X can serve as a rational basis for 
structural design. It is, in a sense, an envelope of ex- 
pected most severe stress in the elastic body. It must 
be remarked that such a stress distribution does not 
occur in any particular one loading but is only the most 
probable value for a large number of loadings. 

From Eq. (108) it is seen that the standard deviation 
of the largest value of stress is also proportional to 
Vv 7 (cf. Eq. (92)]. Hence, following the example in 
Section (VIII), the life expectancy is inversely pro- 


ff 
portional to V 7?. 


(X) CONCLUDING REMARKS 


The preceding analyses are based on simplified models 
and are not immediately applicable to practical prob- 
lems encountered in aircraft designs. The lines of 
thinking, however, are demonstrated, and it is not 
difficult to imagine various avenues of applications. 
For example, gust responses and landing and take-off 
loads during taxiing over rough grounds or water sur- 
faces may be regarded as stationary random processes, 
while the disturbances to the landing impact—devi- 
ations from a mean curve, such as those represented by 
the curves on the right-hand side of Fig. 3—may be 
regarded as a nonstationary random process. T 

It is obvious that the calculation of the mean square 
response can be extended to systems having many 
degrees of freedom. Such calculations will show, in 
particular, the scale effect—i.e., the effect of size of 
the aircraft on the dynamic responses (gusts, tail buf- 
feting in fully developed wing wake, etc.). Housner’s 

t The response to the mean curve can be treated as an ‘‘exact”’ 
The present approach is entirely different from an 
See reference 


problem. 
earlier proposal made by Biot and Bisplinghoff.*! 
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“envelope’’ treats a process that can hardly be regarded 
It is directly applicable only to the 
Nevertheless, 


as stationary. 
particular type of excitation discussed. 
it is indicative of the general form of the result to be ex- 
pected with regard to the spatial distribution of 
stresses. It shows in particular, that the equivalent 
static inertia loading, commonly used in current air- 
craft design, can be misleading and unconservative. 

The application of the theory of extreme values to 
the ensemble averages leads to an envelope of the most 
probable largest stress in a structure. Such an enve- 
lope can be used as follows: When the distribution of 
the largest stress over space is known, it is necessary 
only to measure experimentally the intensity of the 
largest stress at a few points on the structure. Such 
measurements, taken over a long period of time, will 
determine the time interval in which a specified maxi- 
mum intensity of stress will occur once, in association 
with a given probability. Thus the life expectancy of 
the aircraft can be determined for a given specification 
of the largest stress. Conversely, with the service 
life of the aircraft specified, the specification of the de- 
sign load can be made. The reasoning is analogous to 
the determination of the maximum load factor by the 
V-G diagram. 

It is important to realize that the stress envelope re- 
lates to a large number of experiments and not to any 
specific forcing function. If the design stress is based 
on such an envelope, it will ensure, im the long run, a 
uniform factor of safety over the entire structure. 

As was shown above, the solution of the structural ad- 
mittance requires the solution of aeroelastic equations 
similar to those governing flutter, dynamic stability, 
etc. Hence, it is clear that the proximity of the flight 
speed to any of the critical speeds of aeroelastic in- 
stabilities will induce excessive dynamic stresses. 

Some additional results relating to less restrictive 
models are available. However, too brief a presenta- 
tion of quantitative results might prevent an intelligent 
interpretation. Hence, these results will be reserved 
for future publications. 
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An Experimental Study of Natural 
Vibrations of Cantilevered Triangular Plates’ 


P. N. GUSTAFSON,? W. F. STOKEY,? anp C. F. ZOROWSKIi 
Carnegie Institute of Technology 


ABSTRACT 


Experimental results are obtained for the lowest six natural 
frequencies and the associated nodal lines of triangular plates of 
uniform thickness clamped on one edge. Two series of plates are 
examined. In one series, a delta-wing plan form of fixed span is 
maintained, and five plates are examined having aspect ratios 
ranging from 2 to 20. In the other series, a fixed span and a 
fixed root chord are maintained, and five plates are examined 
having mean chord-line sweepback angles ranging from 0 to 45°. 
Test results are presented by graphs and photographs. These 
results show that in the delta-wing plan form of fixed span, the 
frequency of each of the six modes increases as the aspect ratio 
increases No such 


general property is exhibited by the plates in the sweepback series, 


that is, as the root chord is made smaller 


except for the ‘‘bending’’ modes, where the frequencies decrease 
with increasing angle of sweepback. A simple method is de- 
scribed for finding the natural frequencies of a plate of any size, 
thickness, and material as long as its plan form is similar to one 
of the plates used in the experiments, or any other plate whose 


frequencies are known. 


INTRODUCTION 


_ UNIFORM FLAT TRIANGULAR PLATE Serves as an 
approximation to rudimentary wing and tail sur- 
faces on rockets, bombs, and some aircraft. There is 
also some evidence that flat-plate vibrations may be 
correlated with actual delta wings of hollow tapered 
construction. 

The determination of the natural modes and fre- 
quencies of any thin flat plate may be reduced to the 
following boundary-value problem:! 


AAw = (¢/D) »w (1) 
with the boundary conditions 
ow , 
= 0 (2) 
= On 
along a clamped edge, and 
(0°®w/On*) + (2 — pw)d*w/Onds? = 0 
Ow wd'w 0 (3) 
On* Os” 
along a free edge, where 
Received November 5, 1952 


* This report is based on research conducted in part under a 
contract between the Flight Research Laboratory of the U.S 
Air Force and Carnegie Institute of Technology 

+t Department of Mathematics. 

t Department of Mechanical Engineering 


transverse displacement 


aw 

o = mass per unit area 

D = flexural rigidity = Eh*/12(1 — yu?) 

A = square of angular frequency 

n = normal coordinate at any edge (4) 
s = tangential coordinate at any edge 

k = modulus of elasticity 

u = Poisson’s Ratio 

h = plate thickness 


There is no exact analytical solution available for this 
problem at the present time. Some success has been 
achieved with Rayleigh and Rayleigh-Ritz methods, 
but thus far only a few modes and a few plate shapes 
have been examined. There appears to be a great 
need for experimental solutions of this problem, not 
only to supply immediate answers but also to serve as a 
check for approximate analytical solutions as they are 
developed. 

It was decided to explore the effects that plate shape 
might have on modes and frequencies by studying the 
separate effects of the two parameters necessary and 
sufficient to describe the plan form of any triangular 
plate—namely, aspect ratio and sweepback angle. 
The usual definitions are assumed but will be repeated 
for clarity, as shown in Fig. 1. The plate configura- 
tions adopted for the experiment and their designations 
are given in Figs. 2 and 3. 

Note that plates A2 and S3 are identical so that only 
nine different plates are represented in the two groups. 
All of the plates in the ‘‘S”’ series have a fixed span and 
a fixed aspect ratio of 4.0. The sweepback angles vary 
from 0 (isosceles triangle) to 45°. However, the plates 
in the ‘‘A’’ series do not have a fixed sweepback angle. 
Rather, they are all of the delta-wing shape—that is, 
the trailing edge is a straight line from wing tip to wing 
tip. This represents a concession to the widespread 
interest in such configurations. 
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Equipment for measuring frequencies and indicating 
nodal lines 


DESCRIPTION OF TESTS 


Equipment. Fig. 4 is a photograph showing one of 


the plates being tested. The plates were cut from 


sheet steel averaging 0.06] in. in thickness. Each 
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plate had at its base a 6-in. projection for mounting 
The support consisted of a '/4-in. steel plate to which 
the plate projection was clamped by another piece of 
'/4-in. steel plate and a length of 6-in. steel channel 
The entire assembly was held together and clamped 
firmly to a heavy concrete base by means of two |, 
in. studs embedded in the concrete and three ! ,-in 
bolts anchored in expansion nuts in the concrete. 


The method of exciting the natural vibrations is es- 
sentially that described in reference 2. The plates were 
vibrated by a U-shaped electromagnet to which power 
was supplied by a 25-watt amplifier. The amplifier 
was drived by a Hewlett-Packard Model 200a audio 
oscillator. The magnet was mounted on a_ pedestal 
that could be moved about, since the best place at which 
to drive the plate depended on the shape of the plate 
and the mode of vibration. In particular, it was im- 
possible to excite a mode with the magnet directly be- 
neath a nede line for that mode. The magnet mount 
was equipped with a screw mechanism for vertical ad- 
justments. 

A microphone was placed above the plate and its out 
put fed to the Y-input of a cathode-ray oscilloscope. 
This microphone picked up the sound emitted by the 
vibrating plate, and, since the intensity of the sound 
depended on the amplitude of vibration, the oscillo 
scope served as an indicator of changes in amplitude. 
Since the magnet pulled twice on the plate for each 
cycle of the oscillator, the plate usually vibrated at 
twice the oscillator frequency. However, in certain 
cases it was found that the plate could be driven at 
some natural frequency by setting the oscillator at a 
value other than one-half the plate frequency. For 
instance, the fundamental frequencies were excited 
in most cases by setting the oscillator at the natural fre- 
quency. This made it necessary to have some means 
of determining the ratio of plate frequency to oscillator 
frequency. This was done by feeding the oscillator 
output to the X-axis of the oscilloscope to which the 
microphone was connected and observing the Lissajous 
patterns. A further check was made where possible by 
measuring the plate frequency directly with a strobotac. 
Power-line frequency was used as a standard in all 
tests, the oscillator being calibrated against this over 
the range of frequencies used. 

Procedure.—The plate to be tested was clamped 
firmly in the support, and the magnet was placed be- 
neath it. 
the plate vibrated at its fundamental frequency, which 
The oscillator 


The oscillator frequency was increased until 


was indicated by a large amplitude. 
frequency was next increased gradually until resonance 
occurred at the second mode. This again was indicated 
by a marked increase in amplitude and also by an in 
crease in the intensity of the sound emitted, which 
could be heard and was also observed on the oscillo- 
scope. This same procedure was followed for each 


higher mode in turn up to six modes. 
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NATURAL VIBRATIONS 


The closest possible estimate of the natural fre- 
quency was determined in each case by reducing the 
power input to the magnet to the lowest value at which 
vibration could be observed, at the same time carefully 
adjusting the oscillator frequency to give maximum vi 
bration as observed by the microphone output on the 
oscilloscope. This technique narrowed the range of 
frequencies at which resonance could be detected. 
Some of the plates had fundamental frequencies lower 
than the lowest oscillator frequency. For such cases 
the fundamental frequency was determined by pluck 
ing the plate and measuring the frequency with a stro 
botac. 


The nodal lines were determined by the well-known 
method of sprinkling the plate with a light granular ma- 
terial. Ordinary table salt proved satisfactory and af- 
forded sharp contrast with the plates, which were 
painted with a thin coat of black lacquer. When the 
salt was sprinkled over the plate, the crystals either 
went off the edge or settled along a nodal line. These 
lines were photographed after the power input was re 
duced and the vibrations had damped out. The fre- 
quency values presented in this report were determined 
before the plates were painted. It was observed that 
the paint altered the natural frequencies slightly but 
had no noticeable effect on the nodal lines. 
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Second, third, and combined modes of 10- by 10-in 
isosceles plate 


Fic. 5. 


An interesting phenomenon occurs when two succes 
sive natural frequencies of a plate are close together or 
when a natural frequency is nearly an integral multiple 
of another. In either event the plate may vibrate in a 
superposition of the two individual modes,* and the 
only point (or points) permanently at rest are therefore 
at the intersection of the nodal lines associated with the 
individual modes. In fact, it is sometimes difficult 
to make the plate vibrate in the individual modes con- 
cerned without artificially restraining it at several points 
along a nodal line. This behavior is shown in Fig. 5, 
which gives the second mode (166 cycles per sec.), 
third mode (169 cycles per sec.) and combination of the 
second and third modes for the 10- by 10-in. isosceles 
plate. 

In attempting to find the sixth mode of plate A5, 
a series of four node points was observed, indicating a 





Nodal lines and frequencies for delta-wing series. 
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Fic. 7. 


superposition of the fourth mode (608 cycles per sec.) 
and the sixth mode (1,225 cycles per sec.). In this 
instance, the sixth-mode frequency was nearly double 


the fourth-mode frequency. 


DISCUSSION OF RESULTS 


Delta-Wing Series.—Fig. 6 shows the photographs 
of the nodal lines for the first six natural modes of the 
five plates in the delta-wing series, arranged in order 
The frequency of each mode 
Definite trends in both 


of increasing frequency. 
accompanies the photograph. 
frequency and mode shape may be observed if one 
selects a given mode of the widest plate and traces 
through its behavior to the narrowest plate. In each 
case the frequency increases as the chord becomes nar- 
rower. The fundamental mode of all five plates is 
primarily a ‘“‘bending’’ mode, resembling the behavior 
However, some chord- 
For the 


of a simple cantilever beam. 
wise variation in displacement is present. 
wide plate, Al, the fundamental is the only mode of 
the six examined which can be given any such simple 
classification. But as the plates become narrower they 
approach more nearly the behavior of a beam, until in 
the plate designated A5 all of the modes are roughly 
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Nodal lines and frequencies for sweptback-wing series. 


speaking ‘‘bending’’ modes except the fourth, which is 
a ‘‘torsional’’ mode. 

A better picture of the variation of frequency with 
aspect ratio is afforded by the graphs of Fig. 8, where 
frequency is plotted against aspect ratio for each of the 
six modes. If we are to expect a smooth behavior from 
the curves, it is essential to fix our attention on those 
vibration modes belonging to the same “‘shape family,” 
as determined from an inspection of the nodal lines of 
Fig. 6. Family “A” 
mental modes of all five plates. Family ‘“‘B’’ consists 
Family ‘“‘C’’ consists of the 
third modes of plates Al, A2, A3, and A4, and the 
That is, the mode shape that 


on Fig. 8 consists of the funda- 
of all the second modes. 


fourth mode of plate A5. 
starts out as the third mode of the lower aspect ratios 
becomes the fourth mode of the high aspect ratio plates, 
and vice versa, at an aspect ratio of approximately 12.9, 
as seen on Fig. 8. Families ‘“E’’ and “‘F’’ consist of 
the fifth and sixth modes, respectively, for all plates, 
so that in this series there is only one instance of an 
interchange of frequency. 

Sweptback Wing Series.—Photographs of the nodal 
lines for the first six modes of each plate in the swept- 
back series are shown in Fig. 7, with the modes arranged 














FREQUENCY IN CYCLES PER SECOND 


Fic. 8 


in or 
wher 
back 
the | 
As in 
ful in 
belon 
mode 
famil 
wher 
and | 
essen 
ys eh 
in Col 
the fr 
matel 
be m: 
it is 
plate 
may | 
The f 
off as 
ing’ 1 
a fixe 
increa 
conse 
famili 
to inc 
to tan 


Con 
this e 
differe 
case ( 
vibrat 
examy 





h is 


vith 
lere 
the 
rom 
10Se 
ly,” 
s of 
da- 
ists 
the 
the 
hat 
HOS 
tes, 
2.5, 

of 
Les, 

an 


dal 
pt- 


red 


4 











NATURAL VIBRATIONS 
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Frequency vs. aspect ratio for plates of delta-wing shape. 





Fic. 8. 
in order of increasing frequency. The graphs of Fig. 9, 
where frequency is plotted against tangent of the sweep- 
back angle, show clearly how the frequency changes as 
the plates are swept back through increasing angles. 
As in the case of the preceding series, one must be care- 
ful in plotting these curves to consider only mode shapes 
belonging to the same family. With one exception, all 
modes of a given order also belong to the same ‘‘shape 
family.’’ That exception is in the isosceles plate, SI, 
where the third mode is essentially a ‘‘bending’’ mode 
and belongs in family ‘‘B,’’ while the second mode is 
essentially a torsional mode and belongs in family 
“C.” These are the two modes previously discussed 
in connection with Fig. 5. The point of interchange of 
the frequencies occurs at a sweepback angle of approxi- 
Several general observations may 
First, 


mately tan—! 0.035. 
be made from inspection of Fig. 7 and Fig. 9. 
it is again noted that the fundamental mode of each 
plate is primarily a ‘“‘bending’’ mode. Family ‘‘B”’ 
may also be roughly classified as a ‘‘bending’’ mode. 
The frequencies in families “‘A,’’ ‘‘B,’’ and ‘‘D’’ drop 
off as the sweepback angle is increased. 
ing’’ modes, this may be explained by the fact that, for 
a fixed span, an increase in the sweepback angle also 
increases the effective length of each wing half, with a 
consequent drop in frequency. For the other mode 
families ‘‘C,’’ ‘“‘E,”’ and ‘‘F,”’ the frequency is observed 


For ‘‘bend- 


to increase up to sweepback angles of from tan~'! 0.3 
to tan~' 0.5 before dropping off. 


Contrast Between Plates and Beams.—The results of 
this experiment serve to emphasize several important 
differences between beam and plate vibrations. In the 
case of a cantilever beam the natural frequencies of 
vibration are, on the whole, widely separated. For 
example, in a uniform cantilever beam, the first five 
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bending frequencies are in the ratios of 1, 6.27, 17.55, 
34.4, and 56.8.4. On the other hand, the natural fre- 
quencies of flat plates are much closer together, as for 
example in the 10- by 10-in. right triangular plate, where 
the first six frequencies are in the ratios 1, 3.94, 5.50, 
9.42, 12.8, and 16.74. It would appear that approxi- 
mately fwice as many modes of plate vibrations must be 
examined to cover the same frequency range as in beam 
Probably the same observation applies in 
This 


vibrations. 
comparing delta wings with conventional wings. 
suggests that, if three natural modes are required in the 
flutter analysis of a conventional wing, then six natural 
modes would be needed for a comparable flutter analy- 
sis of a delta wing. 

Another point of difference arises in connection with 
the number of node points and nodal lines. In a canti- 
lever beam, each higher bending mode is accompanied 
by the addition of one more node point. However, 
in the case of flat plates there is no corresponding rela- 
tionship involving the nodal lines. In fact, many cases 
exist where a higher mode actually possesses fewer 
nodal lines than a lower mode, as evidenced by the fifth 
and sixth modes of plate Al or the third and fourth 


modes of plate A5. 


THE FREQUENCIES OF PLATES OF SIMILAR PLAN FORM 


In order to make use of the experimental data, it is 
necessary to convert from the frequency of any given 
plate to the frequency of a plate of an exactly similar 
plan form but having different size, thickness, and ma- 


terial properties. Suppose the frequency of plate (A) 
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in Fig. 10 is known, and it is desired to find the fre- 
quency of plate (B) in which all dimensions except 
thickness are some scale factor a times the correspond- 
ing dimension of plate (A). 

In the following development, primed quantities will 
refer to plate (B) and unprimed quantities to plate (A). 
Now suppose the exact displacement, w(x,y) is known 
for a vibration mode of plate (A). This function will 
satisfy the differential Eq. (1) and also the boundary 
conditions (2) and (3). 

Next one might predict that, in the corresponding 
natural mode of plate (B), the displacement would be 
exactly similar to plate (A), and since the magnitude of 
displacement has no effect on frequency (as long as all 
displacements are small), one might guess 


w'(x’,y’) = w(x,y) (7) 
provided that 
x’ = ax | 
> (8) 
y’ = ay 


We then inquire under what conditions the function 
w’ will satisfy the differential equation 


A’A’(w’) = (0'/D’)\X'w’ (9) 


with the boundary conditions 


w’ = Ow’/dn’ = 0 (10) 
along a clamped edge and 
(O*w’/On’?) + (2 — pw’)(O0%8w’/Onds’?) = 0 \ (11) 
(0°w’/On’?) + p’(0?w’/ds’2) = 0 f 
along a free edge. Now 
Ow’ Ow’ dx 1 Ow 
= = - (12a) 
Ox’ Ox dx’ a Ox 
ow’ od (2) dx 1 Ow we 
Ox’? = Ox \Ox’/ dx’ —s a? Ox? _— 
Similarly, 
O'w’ = 1 Ow O'w’ —s 1 O*W 
Oy"* a® Oy* On’? — a On? 
ow’ = 1 Ow Ow 1 dw 
On’ aon On’ds’? a Onds? | ‘i 
r (13) 
O'w’ = 1 O*w AVA" (wy! Ls 
On’? = a On? “= a‘ ™ 
O7*w’ _ 1 Ow 
Os’? a? Os? 


Substituting Eq. (7) into Eq. (9) and using Eqs. (13), 
we obtain 
(1/a‘)AAw = (0'/D’)\'w 


or 


AAw = a‘(o’/D’)\'w 
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By comparison with Eq. (1), this will be satisfied if 


o o 


‘ i, d 
"a? “s 
or 
1D’ oa , 
> A (15) 
at D ao’ 


The boundary condition (10) is also satisfied upon sub- 
stitution of Eqs. (7), (2), and (13) into Eq. (10). Turn- 
ing now to boundary conditions (11) and using the rela- 
tions (13), we obtain 


O*w’ ; fw’ ms 
on’s *e~ 7 On’ds’2 


l & O*w (16 
nt nce (2 = 2) = )) 
a’ Lon? + , Onos? 


44 ow si | w ; | (17) 

“ Os’? a? LOn? Te Os? i. 
It is clear from Eq. (3) that the right-hand sides of 
Thus, 
although a similar displacement will satisfy the bi- 
harmonic equation, it will not satisfy the free-edge 
boundary conditions unless the two plates have equal 


O7w’ 
On” 


Eqs. (16) and (17) will vanish only if up’ = up. 


Poisson ratios. However, most metals do not differ 
by much in this respect, so that Eq. (11) should be 
closely approximated by similar displacement functions 
Then Eq. (15) may be 
used to compare the two frequencies and may be written 


even for different materials. 


as 
1 E’h"® (1 — pw?) ph 
at Eh® (1 — p’?) p/h’ 
1 E’ph”?(1 — p? 
e ph’?(1 pe”) (18 
at Ep’h?(1 — w’?) 
where 
h = plate thickness (19 


p = mass per unit volume 


9 


Or, since A = wo? = 


ie ee a= (20 
’ ah LE p’ (1 — pw”) 


(2xf)?, 
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NATURAL 


To summarize, Eq. (20) may be used to convert fre- 
quency from a given plate of one size to any other plate 
of similar plan form but of different size, thickness, and 
material. The transformation is exact for plates of the 
same Poisson ratio and provides a good approximation 
for plates having Poisson ratios not greatly different. 


CONCLUSION 


Several recent papers®~’ deal with natural vibrations 
of rectangular cantilevered plates, considering the effects 
of aspect ratio and sweepback. The same two param- 
eters are considered in this report, but in the aspect 
ratio series a delta-wing plan form is preserved, which 
results in a variable sweepback as the aspect ratio 
changes. A more systematic (and more comprehen- 
sive) program would consist of the selection of, say, 
five sweepback angles, as in series ‘“‘S.’’ For each angle 
The delta-wing plan form 
Then 


” 


examine five aspect ratios. 
could be included as a special case for each angle. 
25 plates would be required in all, as compared with 
nine covered by this report. There is also the clipped- 
wing plan form to consider, and possibly a simple corre- 
lation may be found which would enable one to predict 
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truncated plate vibrations (trapezoidal plan form) 
from the vibrations of the complete triangular plate. 
These progams may be undertaken if a sufficient de- 
mand should arise. 

The authors hope that this report will prove useful 
to design engineers and also to those interested in com- 
paring analytical solutions with experimental evidence. 
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Some Effects of Surface Curvature on 
Laminar Boundary-Laver Flow 


JAMES S. 


MURPHY? 


Unwersity of Michigan 


SUMMARY 


The laminar flow of a viscous incompressible fluid over a two 
dimensional curved surface is investigated for two cases, one in 
which the curvature is ‘“‘large’’ and the other in which it is 
“‘moderate.’’ The boundary-layer equations applicable to these 
cases are obtained as approximations from the exact equations of 
motion by an order-of-magnitude analysis. These equations are 
solved for flow over a particular surface with zero surface pressure 
gradient. In this analysis, the pressure gradient normal to the 
surface is included, and the outer boundary conditions are modi- 
fied in accordance with the requirements of flow over a curved 
surface. 

The results indicate that for equal Reynolds Numbers, the 
stress on convex surfaces is less than the flat-plate value, while 
The 
most important effect of surface curvature, for the cases con 
sidered, is the modification of the shape of the velocity profile 
The requirement 


the stress on concave surfaces is greater than for a flat plate 


near the ‘‘outer edge’ of the boundary layer. 
that a smooth transition exist between the viscous flow and the 
potential flow at the outer edge of the layer causes the profile to 
have a negative slope near the outer edge for convex surface 
curvature and a positive slope for concave surface curvature. 
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Fic. 1. Coordinate system for two-dimensional motion past an 


arbitrary cylindrical surface 

EQUATIONS OF MOTION 
ONSIDER THE MOTION OF A viscous incompressible 
fluid along a curved two-dimensional surface. Let 
x represent distance measured along the surface from 
the stagnation point and y represent the distance nor- 
mal to the surface, as shown in Fig. |. Designate the 
velocity components in the directions x and y by u and 
v, respectively. Also denote the curvature of the sur- 
face by A = A(x), which is a continuous function of x 
and is positive for convex and negative for concave 
curvature. Then for steady motion the Navier-Stokes 


equations may be written (see reference |) as 
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where Uy) and L are characteristic velocity and length in any flow under consideration. 


dimensional quantities the equations of motion are 
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In terms of these non- 


| , Ou’ , Ou’ K'u'v’ I Op’ l l Ou’ = Ou’ 
r , u , ss v , ~ ai , iis _ re / + € , a — 
i+ K’y ov oy 1+ K’y 1+ K’y’ ox’ RL(1+ K’y’)? Ox”? = Oy” 
y’ OK’ ou’ K’ ou’ K"*u’ 2K’ Ov’ 1 y’ | 4 
(1 + K’y’)® Ox’ Ox’ (1 + K’y’) oy’ (1 + K*y’)* (1 + Ky’)? Ox’ (1 + K’y’)* Ox’ a 
1 Ov’ Ov’ K'u"? Op’ l l O?v’ = O?v”’ 
ag to - a - + saat sn 
1+ K’y Ox oy 1+ K’y oy RL + Ky’)? Ox"? oy”? 
y’ OK’ Ov’ K’ Ov’ ted u’ OK’ 2K’ y’ 
_ (4b) 


(1 + K‘y’)® Ox’ Ox’ 1+ K’y’ ody’ 


where 


R = pU.L/p = 


and the equation of continuity is 


(Ou’/dx’) + (d/dy’)[(1 + K'y’)v’] = 0 


BOUNDARY-LAYER EQUATIONS OF MOTION 


The above exact equations of motion are extremely 
dificult to solve. Indeed, only a few such exact solu- 
tions have been found. Hamel one by re- 
stricting the stream lines of the viscous flow to coincide 
The resulting viscous 


found 


with those of a potential flow. 
motion has stream lines in the form of logarithmic 
spirals.2 Anothor solution of the complete equations of 
motion was obtained by von Karman for the case of 
steady motion about a rotating disc.* However, by 
means of the order of magnitude analysis‘ we may 
simplify the equations of motion and continuity and ob- 
tain the ‘boundary-layer equations.’ The boundary- 
layer thickness parameter 6’ is taken as small, and all 
terms in the equations are arranged according to the 
hierarchy of orders of magnitude. 


i a 


O 
57 <1 512 


... <0(5'2) < O(6') < O11) < O 
(5) 


Physical experience also shows that the following as- 
sumptions are valid for boundary-layer flow, 


OCI) 
O(1) 


Ou’/dx’ = 
Op’ /Ox’ 


‘= O(1), 
O(1), 


Uu 
02u'/dx’? = 


(6) 


where the primed quantities are nondimensional. Then 
moderately curved surfaces may, in the spirit of the 
foregoing definition of order of magnitude, be charac- 


terized by the condition 


K' = 0(1) 


(7a) 


while surfaces with large curvature may be said to be of 
order 


(7b) 


~ G4 ny’)? - 


(1 + K’y’)® Ox’ (1 + K’y’)? Ox’ 


Reynolds Number 


(4c) 


The conditions (7a) and (6) when applied to the con- 
tinuity equation indicate that the orders of magnitude 
and its derivatives are 


, 


of the velocity component v 


; , Ov’ 
y = 0(5'), , = O(1), 
Oy 
Ov’ 0’y’ 
_ = 0(8), ag = ((1) (8) 
Ox Ox ‘Oy 


provided that the derivative of the curvature has an 
order of magnitude no greater than 


0K’/dx’ = O(1/6") (9) 


For example, integration of Eq. (4c) yields 


y’ ra) , 
= [ = dy’ 
J0 Ox Z 


Hence, if the orders of magnitude (6) hold, we may write 
for any point in the boundary layer 


fg O(1) y’ 
v| < oe dy’ = 
1+ K sa i 


Then for a convex surface satisfying the condition (7a), 
(K’ > 0), 


(1+ K’y’)v’ = 


O(6') 
1+ K’y’ 


1 l l 
wry —e" = O(1) 
|\1+ K'’y 1+ |K’y 1+ O(6’) 
while for a concave surface, (A’ < 0), 
1 | I 
j ns ! ns < = O( l ) 
1+ K'y’| ~1—|K’y’] — 1-—06€’) 
so that 
v’ = O(8’) 


In this manner, all terms involving the velocity com- 
ponents and their derivatives which appear in the 
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equations of motion may be assigned an order of magni- 
tude. Inspection of the orders of magnitude of the 
terms in the first component of the equation of motion, 
Eq. (4a), reveals that all terms in the square bracket on 
the right may be neglected in comparison to 0°u’/dy’?, 
which is O(1/6’2). Also the term K'u'v’/(1 + K’y’) on 
the left side may be neglected in comparison to the 
others. All remaining terms are O(1) except the term on 
the right involving the Reynolds Number, which has an 
unknown order of magnitude. Hence, we conclude that 
this term can be of order no greater than unity and find 
that 


1/R = O(6"?) (10) 


This result, applied to the second component of the 
equation of motion, shows that all terms may be neg- 
lected in comparison to the term K’u’?/(1 + K’y’) on 
the left side. Since this term is O(1), we conclude that 


dp’ /dy’ = O(1) (11) 


Hence, the dimensional form of the equations governing 
the motion in the boundary layer on a curved two- 
dimensional surface, which is consistent with the above 
orders of magnitude, is 


P Ou a Ou _— 1 Op ; O*u (12a) 
Ox oy p Ox oy? 
Ku? = (1/p)(Op/oy) (12b) 
and continuity is 
(Ou/Ox) + (Ov/dOy) = 0 (12c) 


when the surface curvature satisfies the conditions, 


K’ = O(1), OK’/d&’ = O(1/58’) (12d) 


and the Reynolds Number, the condition (10). 

In a similar manner, it can be shown that for convex 
surfaces with large curvature, the boundary-layer equa- 
tions of motion are 


l Ou Ou Kuv 
u v = 
1+ Ky Ox oy 1+ Ay 
l Op i 4 K Ou Ku ) 
p(l + Ky) ox oy’ 1+ Ky oy 1+ Ky’ 
(13a) 
Ku?/(1 + Ay) = (1/p)(Op/Ooy) (13b) 
while the continuity equation is 
(Ou/Ox) + (0/Oy)[(1 + Ay)v] = 0 (13¢) 


when the curvature satisfies the conditions 


l OK’ I 
K’'= o( ), — = o( ), ey > 60 (134) 
5’ ox 6’ 


The restriction to convex curvature is imposed by the 
singularity in the equations of motion at Ay = —l. 
As would be expected, the equations of motion for large 
convex curvature are also valid for moderate curvature 
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and should yield, actually, a more accurate solution fo; 
moderate curvature than do the Eqs. (12). Therefore 
Eq. (13) may be applied to both concave and convex 
surfaces when the curvature is moderate. 


BOUNDARY CONDITIONS 


The no-slip condition at the surface requires that the 
velocity vanish at y = 0. This yields the two boundary 


conditions 
u=v=0 (14) 


at y = 0. In addition, a smooth transition is required 
between the motion of the fluid in the boundary layer, 
and the motion represented by the potential flow in the 
region exterior to the layer. Thus, it is necessary that 
at the ‘‘outer edge” of the boundary layer, the solution 
of the boundary layer tends asymptotically to an ex- 
pression representing approximately the potential flow 
in the neighborhood of the surface. This approximate 
representation characterizing the nature of the poten- 
tial flow near the surface depends upon the magnitude 
of the surface curvature and therefore will be different 
for the two cases considered here. In both cases, how- 
ever, it is necessary to consider only thin boundary 
layers. 

It: is supposed, as in the case of the flow within the 
boundary layer, that in a neighborhood of the surface 
the potential flow is such that 


7 


Ou’ Ou 
- = O(1), >, = OU) (15 
x Ox" 


where the nondimensional quantities are introduced as 
before. Then it follows from the equation of continuity 
that, in the neighborhood of a surface having moderate 
or large curvature, in which y’ = O(6’), the potential 
flow satisfies the further order relations, 


a)! , 


Ov Ov 
,= O18), ,= O(1) (16 
ce) 3) 


x yy 


v’ = O(5'), 


We therefore seek an approximate representation of the 
potential flow consistent with these order relationships. 
Since the potential flow is irrotational, the components 
of velocity satisfy a condition that, in terms of the 
curvilinear coordinate system, can be written in the 
nondimensional form 


l vs On’ K'u’ , 

1+ K’y’ ox’ Oy’ l+ K’y’ Ne 

If the boundary layer is extremely thin, Eq. (17) 1s 

valid at y’ = 6’. For the case of moderate curvature, if 

terms of O(6’) are neglected in Eq. (17), the condition 

of irrotationality may be expressed in dimensional form 
as 


(1/u)(Ou/oy) = —KA (18 


This differential equation requires that 
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Family of surfaces considered in particular solution. 


Fic. 2 
u = C(x)e—*" (19) 
where the arbitrary function C(x) is determined by the 
value of « prescribed by the potential flow at the sur- 
face. If, on the surface, wu = u(x) then an approximate 
representation of the potential flow in the neighborhood 
of a surface with moderate curvature is given by 
u = Uo(x)e—*¥ (20) 

Thus for convex curvature (AK > 0), 

0u/Oy),-5 < 0 

while for concave curvature (K < 0), 

Ou/Ov|,-5 > 0 
In the case of large convex curvature, the approxi- 
mate condition of irrotationality expressed in dimen- 


sional form is 

1 Ou K 

u Oy 1+ Ky 
This equation requires that 

u = D(x)/(1 + Ky) 
where the arbitrary function D(x) is determined by the 
value of u prescribed by the potential flow at the surface. 
If, as before, u = uo(x) on the surface, then an approxi- 


mate representation of the potential flow near the sur- 
face is 
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u = uo(x)/(1 + Ky) (22) 


Hence, for the case of moderate surface curvature, the 
outer boundary condition for the solution of the 
boundary-layer equations is taken as 


u—> u(x)e*", y>O (23) 


For the case of large convex surface curvature, the 
outer boundary condition is 


Uo(x) 
of 7 


a 9) 
1+ Ky ° (24) 


u 
Perhaps it should be noted here that these outer bound- 
ary conditions can also be developed by a combination 
of the second component of the equation of motion and 
the Bernoulli equation. 


FLiow ALONG A PARTICULAR SURFACE 


In order to determine the influence of normal pressure 
gradient and the outer boundary condition, developed 
above, on wall stress for laminar flow, we shall consider 
flow along a particular surface. The effects of surface 
pressure gradient will be eliminated by considering an 
idealized flow past the surface, such that the potential 
flow velocity at the surface is a constant, chosen to be 
U;. 

This flow would be difficult to obtain experimentally, 
but a theoretical solution shows the effects of surface 
curvature and allows a comparison to be made with the 
Blasius solution® for flow along a flat surface. 

The family of surfaces we shall consider has a curva- 
ture distribution given by 

K = A(U,/v)'?x™ 


I/o ~ 
(25) 
where A is a constant of proportionality, U is the con- 
stant potential flow velocity at the surface, and » is the 
coefficient of kinematic viscosity. We shall seek a solu- 
tion for those surfaces for which 


—0.05 < A < +0.05 (26) 


Several surfaces of the family for U; = 100 ft. per sec. 
and v = 0.0001567 are shown on Fig. 2. 

Consider Eqs. (13), which are valid for large curva- 
ture, with the boundary conditions (14) and (24). 
Eq. (13c) may be satisfied by introduction of a stream 
function wy, such that 

oy 1 oy 


= dy’ oe (27) 


, 1+ Ky dx 


Now make the following transformation of independent 


variables: 


y (u\'" 
n=i(“), gas (28) 
2 \vx 
We shall seek a solution for which y has the form, 
y = (U,vt)'"f(n) (29) 
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where f(n) is a function of 7 alone. Then 
1 Uf’ Ou 1 U, f : 
“= ; = — | 
 * oe oe 
Ou =U, (“) : ir 
oy 4 \vé 
O7u a U; (/ ‘) ft ( (30) 
oy? 8 \vé/- | 
"91 + 24m) \ E site i 


Ky = 2An ; 
Substitution of these values into Eq. (12b) yields 
Op pA U*f” (31 
= 31) 

On 2(1 + 2An) 


and using this result in Eq. (12a), along with Eq. (30), 
gives the single equation of motion, 


fi’ + ff’ + 2A ff’ + 2Anff”’ -+- 4Anf”’ 
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for the condition 


Op O& = 0 (33 


2] 


which is required for a constant potential flow velocit, 
along the surface. 


The boundary conditions (11) become 


7 = bad = 0, 7 = O (34a 


and the condition (24) at the outer edge of the boundary 
layer is 


») 


| is 1 + 2Ay’ 


n> 0 (34b 


Assume a solution of (32) of the form 


where the C; are constants. In view of the conditions 


(34a), the first two coefficients are zero. Then substi- 


tuting Eq. (35) in Eq. (32) and expanding and equating 


+ 4A°n?f’"’ + 2Af’"’ + 4A%mf’’ — 4A?f’ = 0 (32) coefficients of powers of 7» to zero gives 
Co CG = —24¢6, C, = 12A°C, 
Cy = —(C? + 96A4°C;) Cs = 10AC,? C; = —128A?C,? — 11,520A°*C, 
Cs = 11C,? + 18324*C,? + C, = —324AC,' — Cio = 8,184A?C2* + 531,936A° X 
161,280A °C, 29,520A4C, — C.2 + 46,448,640A42C, (26 
ov) 
2,580,480A !C, 
Cu = —375C.4 — 203,784A°C,3 — 
10,623,744A°C.? — 
928,972,800A °C, 
Hence, the function f has the form 
. 1 ia at = = 2 a . ca (—1)"A” . ~ 1/3 - 
F os C, *Fo(Ce ‘n) oe AG,(C; *n) + C” G2( Ce *n) Sa a + Cj") 3 G,(Ce 7) + — (37) 
where 
(Con)? (Ca'*n)® 11 Ca')® 375 (Ca)! 
Fo Cs y )= = 7. 
iGo 2! sO TT id 
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14! 17! . 
ve 20Ce'")® 10 C2'*n)® ~~ 324(C2')® 22,896 (C2'7"n) 
Gi(C2/*n) = : _ = — 4+ 
(C2) 3! 6! 9! 12! 
2, 854,146(C2''n) © 595,343,S0S8(C2 /’n)!5 (38 
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The series Fy is that obtained by Blasius for flow over a 
flat surface—i.e., a surface for which K = A = 0. 
Therefore the above solution is a double power series in 
» and the curvature that passes over into the Blasius 
solution as A —~ 0. For small values of curvature such 
that A | < 1, it may be expected that the higher 
powers of A will be small and may be neglected. Thus 
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for —0.005 < A < +0.05, we may approximate solution 
(37) by 
(39) 


x = C, *F (Co *n) aa AG\(C2 *n) 


The outer boundary condition specifies the dependence 


of the unknown coefficient C, on the curvature. In 
general, this condition requires that 
C2" Fy’ — ACGy' + A°G:’ — 2.2. + 
(—1)"A’G,,’ 2 
Co —2)/3 _— 1 +4 2An n » 0 (40) 
which can be written in the form 
C,*[ Fy’ — AnE;’ + A%PEs! — 
9) 
... (—-1)"A"n"E,’ +...) > = » n>O (Ala) 
1 + 2An 
where 
E,’ = G,'/(C2”"n)" (41b) 


For small values of A, this condition can be approxi- 
mated by 


C2" [Fo'(Co') — AnEy'(C2"n)] > 


9 


— 42 
1 + 2An (42) 


» ero 


It has been shown!’ that the value of Fy’ for large (C.'/*n) 
is 


Fy/(Co'") = 1.65562, C2'"n > 3.3 (43) 
The approximate value of E,'(C2'n) can be determined 
by summation of the series expansion of this function. 
The value of £,’ computed in this manner is shown on 
Fig. 3, where, for values of (C.'/*n) > 2, the function has 
been approximated by the tangent at (C2) = 2 in 
order to avoid lengthy computation of further terms in 


the series. The equation of the tangent is 


Ey'(Co'"n) = 0.144C,'"9 + 1.002, (C2’'n) > 2 (44) 
Use of Eqs. (44) and (43) in Eq. (42) yields 
C.**[1.65562 — 0.1444 C.''n? — 1.002An] > 

») 
———-, g>3 (45) 
1 + 2An : . 


Evaluation of C. as a function of A has been accom- 
plished numerically for —0.05 < A < +0.05 and 
large n. The results are shown on Fig. 4. Velocity pro- 
files for A = —0.05, —0.02, and +0.05 are shown on 
Fig. 5, along with the Blasius profile for A = 0. 


DISCUSSION OF RESULTS 


Study of the velocity profiles on Fig..& shows that the 
most important effect of surface curvature is the modi- 
fication of the profile shape for large 7. This modifica- 
tion is imposed by the outer boundary condition. The 
profile for convex curvature would appear to be more 
sensitive to an adverse pressure gradient (Op/Ox > 0) 
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than the profiles for concave surfaces. 


points in the layer on a convex surface where 0u/Oy < 0, 


the flow will be considerably less stable to two-dimen- 
sional disturbances than points closer to the surface.® 
The laminar shearing stress at the wall is given by 


T = p(Ou/OyY),—0 (46a) 
which for the particular flow considered here is 
U; Ui 1 2 P 
T= p C2 (46b) 
4 \ vx 


The variation of C, with curvature shown on Fig. 4 in- 
dicates that wall stress for convex curvature is less than 
the flat-plate value, while that for concave curvature is 
greater than the flat-plate value. 

An estimate of the errors involved in Eqs. (13) is 
possible by consideration of the equations of motion 
made nondimensional by introduction of the following 


quantities: 
x POR u 
Q=— yo=R’* =, a = 
. -—* -” £ 
17, U p : . 
7 = R°—, p=——, Ai = KL (47) 
l 0 pl o~ 


where Uy and L are characteristic velocity and length 
and R = Reynolds Number. 

It can be shown by means of this transformation that 
the errors in Eqs. (13a) and (13b) are of order R~! and 
(0K,/0x)R~**, while the equation of continuity is 
exact. Previous solutions in which all terms involving 
the surface curvature were neglected contained, in ad- 
dition to the above errors, errors of order K,R~ ”* (see 
reference 1, page 121). Perhaps the greatest improve- 
ment introduced here, however, is associated with the 
new outer boundary condition that is valid to order Aj, 
whereas previous solutions contained an error in this 
condition of order A. 

A solution for flow along the surface considered above, 
using the more approximate equations of motion, Eq. 
(12), and the corresponding outer boundary condition, 
Eq. (23), has also been obtained. The equation of mo- 
tion is simpler, in this case, being 


f'' + ff’ + 2Anf”? = 0 (48) 
with the boundary conditions 
f= f' =90, = 0 
Js ’ t (49) 
f’ +2—-—4An, > Of 


The solution of Eq. (48) agrees closely with the above 
solution, Eq. (37), for convex curvature, but the results 
for concave curvature indicate a smaller value of wall 
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Fic. 5. Typical velocity profiles for flow along convex, flat, 
and concave surfaces of the family considered in particular solu- 
tion. ' 


stress than Eq. (37). The difference between the two so- 
lutions may be attributed, essentially, to the differences 


between the outer boundary conditions (23) and (24). 


It should also be noted that the numerical values ob- 
tained in the solution presented here could be improved 
upon by the numerical integration of Eq. (32) with the 
boundary conditions (34). 
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Shock-Wave Effects on the Laminar Skin 
Friction of an Insulated Flat Plate at 
Hypersonic Speeds’ 


TING-YI LI? ano H. T. NAGAMATSUt? 
California Institute of Technology 


ABSTRACT 


An approximate theory on the phenomena of interaction be- 
tween the shock wave and the laminar boundary layer on an in- 
sulated flat plate at hypersonic speeds has been formulated. Re- 
sults on the rate of growth of the boundary-layer thickness and 
the rate of decay of the shock-wave strength have been found 
that hold for M,?° + © R,-! = 0(1). A new set of formulas 
for the average skin-friction coefficient, Cr, over an insulated 
flat plate at hypersonic speeds has been obtained. Calculations 
on the basis of the new Cy formulas yield the data shown in Figs. 
5and 6. Contrary to the conventional theory that predicts a 
steady decrease in Cr as M, increases, the present results indicate 


that Cry may increase with 4, at hypersonic Mach Num- 
bers 
NOMENCLATURE 
L = chord of the flat plate 
X% = interaction distance 
Cy = average skin friction coefficient 
M, = free-stream Mach Number 
x,y = Cartesian coordinates (Fig. 1) 
u,v? = velocity components in x, y directions, respectively 
p = pressure 
p = density 
T = absolute temperature 
uw = coefficient of viscosity 
k = coefficient of conductivity 
R = gas constant 
C, = specific heat at constant pressure 
» = constant index of the power law of viscosity 
Ri, = pilix/m, =.free-stream Reynolds Number based on the 
distance x from the leading edge 
Ri = pil /m, = free-stream Reynolds Number based on the 
chord of the flat plate 
6 «= boundary-layer thickness 
Pr = wC,/k = Prandtl Number 
Tw = pw(Ou/OyY)» = laminar skin friction of the flat plate 


Subscripts 1 and 2 refer to the conditions pertaining to the free 
stream and the outer boundary of the viscous layer, respectively. 
to the conditions on the flat-plate sur 


Subscript w refers 


face. 
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(I) INTRODUCTION 


ly RECENT YEARS, a considerable amount of informa- 
tion concerning the skin friction of a two-dimensional 
thin flat plate moving edgewise through the air at high 
speeds has been obtained through the study of com- 
pressible boundary-layer equations.” The usual basic 
assumption is that viscous effects are confined to a 
narrow region close to the plate surface. The flow field 
outside this narrow region is regarded as nonviscous 
and not effected by the presence of the boundary layer. 
However, when the plate is moving at supersonic speeds, 
disturbances due to the viscous effects in the boundary 
layer can drastically influence the outside flow field.* 
Thus, the retarded flow in the boundary layer causes 
the stream lines to turn slightly away from the plate 
surface, and this results in a shock wave extending 
from the plate into the free stream. This shock 
wave, moreover, decays in its strength as it intersects 
expansion from the 
There is, therefore, a close relation- 


successive waves originating 
boundary layer. 
ship between the rate of growth of the boundary-layer 
thickness in the viscous flow region and the rate of 
decay of the shock-wave strength in the outside flow 
field. 
ary layer and the shock wave are so far apart that their 


In the lower supersonic speed range, the bound- 


interaction effects can probably be neglected as in the 
The situation is quite different 
One 


existing theories.” 
when the hypersonic speed range is approached. 
may define the hypersonic flow régime such that 
M)(v2/u2) = O(1) (reference 4). By the familiar argu 
ment, the boundary-layer thickness must be such that 


6/x = O(1 M,) 

Therefore, in the hypersonic speed range, the shock 
wave, which has a wave angle of the order of the quan- 
tity 1/14, is so close to the plate surface that the entire 
region between the shock wave and the plate surface 
should be considered as a viscous flow layer.** Conse 
quently, for the determination of the shock-wave char 
acteristics, viscous effects must be considered,°® and for 
the estimation of the plate skin friction, due care must 
be exercised to account for the shock-wave effects. In 
** Lees" postulates the region between shock wave and the 


plate surface as consisting of hot and cold layers; thus, his 


formulation is somewhat different from ours 
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FIG.1 — THE COORDINATE SYSTEM AND THE FLAT PLATE 





the present paper, an approximate theory of the shock- 
wave effects on the hypersonic laminar friction on an 
insulated flat plate will be presented. 


(II) Bastc EQUATIONS AND ASSUMPTIONS 


The thickness of the shock wave® 7 is assumed in- 
finitesimally thin.t The flow variables in the region 
bounded by the shock wave and the plate surface 
satisfy the compressible boundary-layer equations, 


Ou + Ou dps + re) ( on) (1) 
1 v-=- 
ae ae dx ' dy\"ay 
+ as 0 (2) 
u - pv = y. 
ox * Oy? 


: (C,T) + : Cf 
1 3 v . = 
= oy er 


dpe fe) (: a 9) 

"i ae ae = 
The order of magnitude arguments responsible for the 
derivation of the above equations can be carried through 
in hypersonic flow problems (cf. Appendix of reference 
However, these arguments can 
In the present 


(3) 


1; also reference 8). 
fail near the base of the shock wave.® 
analysis, the shock wave is assumed to start from the 
plate leading edge (Fig. 1). Therefore, the above equa- 
tions do not hold in the vicinity of the origin, 0 < 
x < e. The flat plate possesses, in a certain sense, 
a fictitious curvature because of the presence of the 
viscous layer. It is because of the fictitious curvature 
effects that the terms containing dp. /dx are retained in 
Eqs. (1) and (3). 

The problem is then to seek the solution of the bound- 
ary-layer equations, Eqs.(1)—(3), with the appropriate 
boundary conditions. The conditions of no slip and 
heat insulation are assumed valid at the surface of the 


7 The ratio of the shock-wave thickness Y and the boundary- 

layer thickness 6 can be expressed as® 

X 87 u*® (M,,* + 1) 

6 y¥ +1 p*un*é (M,,* — 1) 
where the superscript * denotes the quantities at conditions per- 
taining to the critical speed of sound and the subscript » denotes 
the normal component to the shock-wave front. For the present 
problem, calculations show that, except in the region of the lead- 
ing edge of the plate, Y/5< 1 if 

M,,* — 1 >> M270 — » 
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flat plate. At the outer boundary of the viscous layer, 
besides the usual conditions, the pressure jump condition 


across an oblique shock wave must be satisfied. Thus, 
one has the following conditions: 
Aty =0,u =v = 0, 
(O7'/Oy)» = 0 (4) 
At y = 6,u = t&, 
(Ou/Oy)s = (O7u/Oy?)s = 0 ) 
poy (ar = y-1 (5 
a vyti\ ‘dé a 
It may happen that the rate of decay of the shock-wave 
strength is such that at x = x» the shock wave already 


degenerates into a Mach wave. When this happens, 
a new phase in the interaction between the outside flow 
and the boundary layer begins so that the entire prob- 
lem must be formulated differently. The present 
problem is therefore to find the solution satisfying Eqs. 
(1)—(5) in the domain 0 < y < 6,€< x < X%. 

In the present analysis, the following simplifying 
(1) The perfect gas law p = 
(2) the air viscosity coefficient 


assumptions are made: 
RpT is assumed valid; 
is a function of temperature only, and the simple law 
pal is assumed; (3) the specific heat at constant pres- 
sure is assumed constant, and the Prandtl Number is 


assumed unity. Then, the particular energy integral 


C,T + (1/2)u? = constant (6) 


which satisfies the conditions in Eq. (4), reduces Eq. (3) 
identically to Eq. (1). Integrating Eqs. (1) and (2) 
across the boundary layer and combining the results 
yield the familiar momentum integral relation, which 
can be written as 


Tip d pots” > pu u 
_ = z f 1 — dy + 
pitty” dx pit” Jo  pette Us 
d D2 é Uu = 
m (! ) f (2 - ) dy (4) 
pit,” dx \py 0 polls 
By Eq. (6), one can show that 


T y¥- 1 es (“)'( ") 
= ] M;2 1 — 
T> + 2 Te uy Uy" 


which for 14; >> 1 can be approximated in the layer 


0< y< 6[l — (A,/5)] 


: at = l Me es (““) ( 7 “) - 
T> 2 T»2 i Uo~ 


Therefore, Eq. (7) can be rewritten as 


as 


+ One may assume that for x > xp the boundary-layer growth 
is so slow that the nonlinear interactions between the shock and 
the boundary layer appear only as higher order effects and are 
xo) can therefore be designated as the ‘‘interaction 


negligible. 
> xo, the ordinary boundary-layer theory 1s 


distance.”” For x 


assumed valid. 
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u 
2 d D> *6[1 — (A:/8)] ‘le 
=. «= ; pr / als 
pil’ (y — 1)M,?2 dx | pi as 7 
Uo 


[1 — (A1/8)] 2 
I (y- 


fans 
ial 5{1 (Ai/6) | = 


76 
fa : 


For MW, >> 1, 41/6 << 1, Eq. (7a) can finally be ap- 


proximated as follows: 


Tx 2 d k f (u/ U2) | 
= dy} — 
piuy" (y — 1) Mi? dx Lp Jo 1 + (u/ue) ~ 


6 d (”) 
(9) 
yM)? dx \pi 


Eq. (9) was first used by Shen® in his research on the 
viscous effects on the shock-wave shape over an insu- 
lated wedge at hypersonic speeds. 


(III) NONLINEAR INTERACTIONS BETWEEN THE SHOCK 
WAVE AND THE BOUNDARY LAYER 


The momentum integral equation, Eq. (9), can be 
treated by the approximation method of Pohlhau- 
sen." '! We shall assume a velocity distribution 


function as follows: 


, A(x) ; 

f(y*) = 2y* — 2y*3 + y*4 + y*(1 — y*)® (10) 
‘ : ; 6° : 

where f = u/u. and y* = y/6. Except for the pressure 


jump condition in Eq. (5), f satisfies the following given 
boundary conditions: 


fl) =1, f') =f’) =0, f(0) =0 


f"(0) = —A = (67/pytt,) (dp2/dx) tT (11) 


In reference 5, Shen used a velocity profile linear in 
y*, thus satisfying only two of the conditions in Eq. 


+ u; and wm are related across the shock wave by the relation 
Us vo dé 1 
={1+ 
My Us dx 
For hypersonic flow, one can take 


tlo/u,; = 1 0(1/M,2) = 1 


LAMINAR SKIN FRICTION 

y) d 2 

n£(A)x 

pil,” dx Pr 

y T> (“) Uo tly + d ¥ 
—l1ldy 

1)M,? T, \ue ' («) : dx 


pitty” dx 


pi 


u (1 = “) 
Uo Ue pi d (*) 
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(7a) 


(11)—viz., f(1) = 1,f(0) = 0. Eggs. (9), (10), and (11) 


can be combined to yield 
d (2 s) __2 (a) dA (? i) ™ 
dx \py 1 — F(A) dx \~i 
w — | DA 
Mm: + FE (2 - ‘) + 78 | 
pil; 1 — F(A) 6 


*1 
F(A) = | dy*/(1 + f) 
0 


where 


and F’(A) = dF dA. Integration of Eq. (12) y 


ps = = _f [1 — F(A)]? X 
pi — F(A)]? Jo 


DA 
(y — | 2—- ) 
’ ( 6 


1 — F(A) 


+ yA | dx 


(12) 


(13) 


ields 


(14) 


The quantity p: p; must, moreover, satisfy the pres 


sure jump condition mentioned in Eq. (5). 
fore, we have the following equation: 


Y (*) -72 yh 
2(y + 1) \dx y¥+1 1M, 


Mw, pilty . > 
[1 - F(A)]? | {1 = F(A) | x 


j(y — 1) [2 — (5A/6)] 1 
{ 1 — F(A) 


where 


l.. 
vA can 


There- 


( 1.) 


(16) 
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On the other hand, the pressure jump condition may be 
substituted in the last condition in Eq. (11) to yield 


another equation 


, y M,*p, dn E l am’ _ 
ae yY+1 Met dx Ldx? 2n & | —" 
The problem is to find A and 7 satisfying Eqs. (15) and 
(17) simultaneously. 7 is the boundary-layer thick- 
ness function, while A is essentially the shock-wave 
pressure function. Eqs. (15) and (17) show the inter- 
dependence of the growth of the boundary layer and 
the decay of the shock wave. The nonlinear interac- 
tion characteristics are clearly exhibited. 


(IV) APPROXIMATE SOLUTIONS OF THE BASIC EQUATIONS 


To deal with the dual relations in Eqs. (15) and (17) 
exactly as they stand, one may eliminate A between 
them to obtain one equation defining 7. The resulting 
equation is, however, nonlinear, a rigorous treatment 
of which would be extremely tedious and difficult. 


2 Dm?’ 7 
n= [Ot | Be | I n— rar (f iy = 
be pill; 0 0 


It can be verified that A = Apo = constant satisfies Eqs. 
(17) and (18) simultaneously. The constant Ag must 
satisfy the following equation: 


Ao = 19) 
in 1 — F(Ao) ' 
Eq. (i8) can be then evaluated as follows: 
| who] 
m=. E +1)° d (20) 
o pil, 


By Eq. (8), #e/: can be expressed as 


w — ] ‘“ 
Me (? us) 
Mi 2 


Thus, the boundary-layer growth follows the following 
approximate law: 


bo 2 1 — l wy 1 2 1 
p. = 7 (¥ + 1) (2 3 ) Ao *M,”* Riz ; (21) 
Ws 


x 


It can easily be shown that 


i, V3 | —1\"* , = 
: = : (y oe 1) , (: ) Ao Me” Bas . 
dx 2 2 


(22) 


In Eq. (22), déo/dx is the approximate slope of the shock 
wave, which, as assumed, forms the outer boundary of 
the viscous layer. The result indicates that as x in- 
creases the slope of the shock wave decreases according 


AERONAUTICAL 


SCIENCES—MAY, 1953 
For hypersonic flow, 1/; > 1, we shall consider the 
following approximate equations instead : 


Y (%") es Mi/ pity x 
2(y + 1) \dx [1 — F(A)]? 


‘ vaxy §(¥ — 1) [2 — (5A/6)] 
1 — F(A)]? < AS de 
J [ ( ] ] 1 _ F(A) + & { dx 


(15a 


i y M,*p, dn EF 1 /dn\? - 
a y+ 1 wet dx Ldx? 2n (77) | oe 
The omission of the term —[(y — 1)/(y + 1)] (nM 
can be justified @ posteriori in a certain limited range of 
x’s [ef. Section (V)]. Notice that some of the nonlin- 
earity properties of Eqs. (15) and (17) have been pre- 
served; the treatment of the approximate system 
should provide an initial step toward the understanding 
of the exact system. From Eq. (15a), 7 is easily ob- 
tained as follows: 


vavye $(Y — 1) [2 — (5A/6)] \ 7 

F(A) }? - Ae dxe Ix 18 
yy 1 — F(A) Ye oe 

tox“. At x = x9, the shock wave will eventually 


be parallel to the local Mach wave direction such that 


(*) es l es V3 + 1 WP ? = ‘) 4 
nie” we” gy YT * 2 


1 2 ltt 
Ao *M,”* & (2) (23) 
Xo 


The “interaction distance”’ x9 can, therefore, be given as 
follows: 


Xo g _ = ‘) ui "s 
= l Ao 
Lit» ( 2 R, 


The approximate law of the decay of the shock wave 


(23a 


can also be expressed in terms of p/p; as a function 
By Eq. (11), it is seen that, for x < xo, 


") : M, 2(w + 1) 
Y-\o R = 


constant 


of x. 
(2) ae pi) ae ae ? —_ 
L/ d(x/L) 4 


Moreover, for the present approximation, Eq. (12) be- 


(24) 


comes simply 


d (BE) (P=) yay MEE 
a 2yAo (25 
d(x/L) \p, L? D) R, 


Combining Eqs. (24) and (25) yields 


ps _ 3Y (? — ‘)’ , x i 
Pp Ay+1)”% 2 ~~ ~~ RY (x/L)' 


(26) 
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Therefore, the shock-wave pressure jump is varying plies that 


1 . - 
inversely to x ’*. As the distance from the plate lead- uM? < [( Ag)? 
SF x ivy" SS Me Pild)).\0 
ing edge increases, the shock wave must become weaker. ™ 
Theoretically, when x = Xo, the shock-wave angle be- In other words, the approximate solutions are valid 


comes the Mach-wave angle and the pressure ratio in such range of x’s that 
»,/p: should become unity. Substitution of %o/Z from A coe ” 
Bq. (23a) in Eq. (26) gives, however, the result Lem Ao/ Ri (<¢) 
Hence, as x — & [cef. Eq. (23a) ], the approximate solu- 
() _ 2y 1.17 tion will be in error. However, results in the last sec- 
Pil 2. 

This conclusion is not surprising because our approxi- 
mate equation, Eq. (15a), is strictly compatible only 

with the approximate shock-wave condition that 


*? 1 , tion do show that the error involved will not be serious. 
An indication of this fact is that the value of p./p; at 
Xo on the basis of the approximate results deviates only 
about 17 per cent from its theoretically true value. 


bs Dy ds\? To provide a remedy, we must consider the exact 
= (ar *) system of equations, Eqs. (15) and (17). As a pair of 

Pi rv ex corrections to the approximate results, we may intro- 
duce the functions Ay and AA so that they satisfy the 


’) CORRECTION FUNCTIONS BY LINEARIZED EQUATIONS , i 
We 2 following equations: 


The approximate solutions in the preceding section . 
; ae Pe eee 2 n = no + An (28) 
are based on the approximate system of equations, Eqs. ek (29) 
(15a) and (17). It is of interest to show in what range 4 — P ‘eit 
of values of x the approximations are valid. This, y dn\? y—1 1m Lae | 
indeed, can be demonstrated by the following consider- Hy+1)\dx) y+t1M? pit; [1 — F(A)]? x 
a By the approximate results of the last section, + , F(A) }2 j(y — 1) [2 — (5A/6)] 4 \| d 
we have — F(A)}? ; : Ag dx 
0 1 — F(A) ba 
yj — | 0 My - ‘ 
4 . : sail | ( ay) : | (30) 
Y iM y* pil / ey” ° ‘ , 
, ¥ M,*p; dn E l (7) | 3) 
"2 A = - —s (31) 
Hu Pitt I [1 — F(Ay) 2 X y+ 1 wet, dx Ldx? 2n \dx 
1 — F(Ao) ]* Jo : - , 
: . Assume that the correction functions are small so that 
{(y — 1) [2 — (5Ao/6)] \ Mu 
+ yAoe dx = 0 Aox the second and higher degree terms in An and AA and 
{ 1 — F(Ao) 7 u . 
a ath their derivatives can be neglected. Thus, we arrive at 


Therefore, the approximate equation, Eq. (15a), im- a pair of dual relations linear in An and AA as follows: 


’ (2) sien y—-1% — Me (20|° ~ OF = Calet y vy faa) + 
ig = C cone 5 


y+1\dx/ dx yt1M2 pw 1 — F(Ao) 
j (a + B) (vy — 1) [2 — (5Ao 6) ]) f ‘ - 
(1 — ZaAo AAdx (32) 
V7 saaedlid 1 — F(Ao) f 0 
/ M,? In) d?*( An) dno 3. (dno\* | d(An) 1 /dnv\?* 
AA = - Y Pi Jno ¢ ~- +| : _ ( "| +4 ( ) (An) 33) 
y+1 pou, (dx dx? dx? 2m \dx dx 2no? \dx 
where a 2y — 1) pr [Cy # 1) Ao(Me/ pitts) | 
YH 1 pylts [7 — (8 + 5a)Ao] 
(dF /dA) 5 (2K 
7 ¢ ¢ J an 4 (34) 30) 
1 — F(Ao) 6[2 — (5Ao/6) ] 
Hence, we have 
Eqs. (32) and (33) can readily be solved. The results , \< V7 
: 7 = My 
are = x ( “+ 1) ) Ay 
. (3 | 2 R, 
S — 1)x° 3 s(y — 1) z 
dy = 5 ls (35) ( ) + ———1 ; 37) 
37(7 — (8B + 5a)Ao|A? i, 3y (7 — (B + 5a)Ao]M;? 
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Ay -— Wf . 8(y — 1) 
A = Ao + - [ -— f+ 5a) o| ae 7 = — i - (49 
7M; 3y(7 _ (Bp + 5a) Ao] M,? 4 
Ao R, x aan i m 
ot ee ' eT (38) are introduced. The slope of the shock wave can they 
- M,? be given as follows: 


— , 2(wt2 
Evidently, for x/L < M,"‘°*”’Ao/Ri, the second term 
in each of the above equations, Eqs. (37) and (38), repre- 
sents a small correction as expected. The importance 


of the correction terms gradually increases as 
x/L > MP?* Ao/R, 
For 
a MPC Tr? 4 


> "=9Q (*) 
b R, ‘ 


the correction terms actually dominate; 
for such cases, the above linearization process should 
fail. Nevertheless, it is interesting to observe that, 
even as x > ©, we get, from Eq. (37), 6/x = 0(1/14), 
which is consistent with the definition of the hypersonic 


to be sure, 


shock régime [cf. Section (I) ]. 


(VI) VALUE OF Ao 


Before we carry on the computations any further, 
we must first determine the value of Ao. In order to 
do this, the function F(A) defined in Eq. (13) must be 
computed for different A’s. Evaluation of the integral 
F(A) has been performed analyticallyt with the results 
shown in Table 1. These values of F(A) have also been 
independently checked by numerical integration meth- 
ods. On the basis of these values of F(A), the value of 


Ao satisfying Eq. (19) is found to be 


TABLE | 
Values of F(A) 


A F(A) A F(A) 
0.2 0.6106 0.8 0.6085 
0.4 0.6101 1.00 0.6087 
0.6 0.6092 


(VII) GROWTH OF THE BOUNDARY LAYER AND DECAY 
OF THE SHOCK WAVE 


On the basis of the result in Eq. (37), the growth of 
the viscous layer downstream follows the following 
law: 


(40) 


6 — (5) [1 4 = (2) ] 
: ” L D\L 


where, for brevity, the symbols 


> 4 lees 
>= ‘lat , 


t Analytical formula for F(A) is given in reference 1. 


jit 
dé ae he : 


dx ~ 4 (*) [: , fay! (43 
L p> 


On taking o = 0, Eqs. (40) and (43) become Eqs. (21 


w4ia 


_~ 
os | & 
etl 
| 


and (22), respectively. For 
0 < (0/3) (x/L)" < 1 


Eqs. (40) and (43) can be expanded as the following 
series: 


we 


dé = > A [ +2e(*) (ox 4 
dx 4(x/L)" 6X\L ari: 
(43a 


It is interesting to point out that in Eq. (47) of reference 
3, Shen obtained a series solution rather characteristi- 
cally similar to Eq. (40a). However, Shen was not able 
to sum his series and obtain a formula like Eq. (40). 
We also see from Eq. (48a) that Eq. (22) underestimates 
the slope of the shock wave. A new “interaction dis- 
tance,’’ Xo, can be found by setting dé/dx in Eq. (48 
equal to 1/1. 

Xo Sl 9 ‘ 2 ‘ 9\! 9 
7? S 2?M,* [2 — 30M,? + (4 — 30M,?)"J]-? (44 


When o = 0, we have x»/L = %/L as given in Eq 
(23a). It can be shown that the pressure jump func- 
tion p2/pi, obtained from substituting Eq. (43) into 
the oblique shock relation in Eq. (5), reduces approxi- 
mately to the form of Eq. (26). To compute (p2/pi)z 
we have the simple relation 


ps = : (2 aia 30M? + (4 = 30M,?)' *) (45 
Pi Zo 2(¥ + 1) 


Notice that 30.1," is, by Eq. (42), a function independ: 
ent of M,. For Ay = 0.91, y = 1.4, the numerical 
value of (p2/p1);, computed on the basis of Eq. (45) 1S 
1.0147. Since this value deviates by only 1.5 per cent 
from the theoretical value of unity, the correction pro- 
cedure devised in Section (V) can be considered as satis- 


factory. Thus, at x = xo defined in Eq. (44), the pres 
ent method yields solutions that fulfill the simul- 
taneous requirements that (dé/dx),, = 1/M/;, and 


(po/ Pr), = 1.0. 
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(VII) AVERAGE SKIN-FRICTION COEFFICIENT ON AN 
INSULATED FLAT PLATE 


In Figs. 2 and 3, curves of %o/L and xo/L have been 
plotted, according to Eqs. (23a) and (44), respectively, 
against M, for R, = 10%, 10’, and 10%. In computing 
these values we take y = 1.4, w = 0.768, and Ap = 0.91. 
It is seen from these figures that the “interaction dis- 
tance’ at the lower Mach Numbers are of insignificant 
proportion; however, as the Mach Number increases, 
the “interaction distance’’ grows rapidly. The rate of 
growth of the ‘‘interaction distance’’ depends on the 
Reynolds Number of the flat plate; at the lower Reyn- 
olds Number, this rate grows noticeably faster. There- 
fore, at a fixed Reynolds Number, the shock-wave effects 
are extended farther downstream as the Mach Number 
increases; at extremely high Mach Numbers, the shock 
wave, so to speak, wraps itself on the entire flat plate, 
x/L > 1. On the other hand, at a fixed Mach Num- 
ber, the viscous effects on the shock wave are much 
more pronounced, as expected, when the Reynolds 
Number is low. To ensure the validity of the initial 
assumption of continuum flow, it is important to impose 
a certain restriction on the relative magnitudes of the 
Reynolds Number and the Mach Number such as the 
criterion of nonslip flow proposed by Tsien'*—viz., 
M,/Ri'* < 1/100. 

The hypersonic viscous flow over a flat plate, accord- 
ing to the present theory, can then be pictured as in 
Fig. 4. Thus, the average skin-friction coefficient on 
the flat plate can be computed by the following formula: 


Cr = (xo/L) Cr, + [1 — (x0/L)] Cr, (46) 


1 Xo 
Cr, = f C,, dx (47) 
Xo J0 
1 71. 
Cr, = / Ce. de (48) 
L— xo Jx, 


where C,, and C;,, are, respectively, the local skin-fric- 
tion coefficients in the hypersonic shock-flow region 
(x < xo) and the conventional boundary-layer flow re- 
gion (xy <x <L).¢ From the present results, we have 


= “2 he o |x 
Cy = ? M :) E (: ) x 
, 2 RET 6U TSA 
w!/2 Xx site o x ~“ 
|» (*) (1 TS y') | (49) 
By von Karman’s results, '* we take 


a a ce: wie ee |) ke 
Ge = 5 VIO a ( ) M, (*) 


(50) 


where 


where V 8f(0) = 1.57. Substitution of these values in 
Eqs. (46)-(48) yields, after simple calculations, the 
following formula: 


t For convenience sake, it is assumed here that Cys of the pres- 
ent theory applies also in the leading-edge region. 
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FIG.2— INTERACTION DISTANCE RATIO vs MACH NUMBER 
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- oo l ‘s ne Ao Xo o Xo 7 a a Ao 
Cr = 2 ( : M,’ “R, z , ag Me + 4- E s (: + ; > 4 
2 bam! 5 a 6 
j ° o od o Xo 0 a Xo ea l ee ent \/s Xo 
loz [ +? Vs ive eae 57 ( . ) M, R, (1 _ \ =) (469 


(v5 zz 
which, for o = 0, becomes 
(46b 


Without the second term on the right-hand side, Eq. (46b) simply gives von Karman’s approximate Cy formula, 


Therefore, the shock-wave effects on Cr are to contribute a certain correction function depending on Ao, M, 
In the above calculations, it is assumed that the hypersonic shock wave is of 


and R, as in Eqs. (46a) and (46b). 
such a strength that it becomes a Mach wave before arriving at the plate trailing edge. 
and (46b) apply, respectively, for the cases xo/L < land %/L < 1. 
Consider, for example, the diagram = xo/ L for the case of R; = 107 in Fig. 3 
Eventually, at 1, = , Wwe find x)/L = 1. The condition of continuum flow is violated when 
In the range of hypersonic eal Numbers, 21 < JJ, < 32, we therefore could anticipate that xy/L > 
For such cases, 


Therefore, Eqs. (46a 


As AJ, increases beyond 10, xp L 
grows rapidly. 
becomes 32. 
1—1.e., the hypersonic shock tends to wrap itself on the entire length of the flat plate (Fig. 4a). 
the average skin-friction coefficient on the flat plate can be simply given as 


1 ie 5 
I | Ce dx (5 


19) in Eq. (51) and reducing yield the following formula: 
1/s Ag o o o a 
(3 (1 +f) +(447 ‘\(Z) x 
» -~ 18) ~_ 
) Ji+2 “( Sie allt +2) |) (51 
\10¢ a ola 
ae Vr\Vz ° N v/]f 


Cr = 


Substituting Eq. ( 


= w 
Cr = 2(? F us) a2 


which, for ¢ = 0, becomes 








; Ao\ (vy — 1\*4 : ay 
Cr = 4V3 (2 =i we) ( = ) (y +1)°“ A, “RR,” uM" (51b 
) ~ 
Eqs. (51a) and (51b) apply, respectively, for the cases vo/L > land x)/L > 1. 
(IX) Discussions ne ‘ y - 
er MACH 
, ; y WAVE 
Before we discuss the numerical results, it seems _ — 
worth while to retrace the development of the present sHoc’ Fs 
° ° ° ° ° ° aaa Vv 
theory in order to bring out its limitations. The hy- Cfs y_ Cty _— > X 
pothesis has been made that the space between the flat- x | 
: se t ° 
plate surface and the shock wave is filled up with a 
=< a > 


viscous laminar flow, which is adequately described by 


the compressible boundary-layer equations. Assum- 


ing, in this viscous layer, a plausible velocity distribu- 
tion function that satisfies the proper boundary condi- 
tions and making use of the pressure jump relation at 
the oblique shock-wave front, we derive a pair of dual 
relations characterizing the nonlinear interactions 
between the shock wave and the boundary layer. A 
pair of approximate relations is then derived in which 
Solu- 
tions are accordingly obtained on the basis of the ap- 
For a fair evaluation of the signifi- 
cance of the solutions thus found, the following ques- 
Are the compressible 


some important nonlinearities are preserved. 
proximate system. 


tions must be answered: (1) 


FIG.4-HYPERSONIC SHOCK PARTIALLY COVERS 
THE FLAT PLATE 
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FIG.4a—HYPERSONIC SHOCK WRAPS ITSELF ON 
THE ENTIRE PLATE 
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FIG.5 — AVERAGE SKIN FRICTION COEFFICIENT ON AN INSULATED FLAT PLATE AT 
HYPERSONIC MACH NUMBERS 


boundary-layer equations valid in the hypersonic viscous 
régime as postulated above? (2) How are the solu- 
tions affected by the choice of the velocity distribution 
function ? 

We have shown that the postulated hypersonic vis- 
cous régime can be defined such that 6/x = O(1/4/)). 
It can be shown (Appendix of reference 1), by the famil- 
iar order of magnitude arguments, that Eqs. (1)—(3) 
are consistent with this definition of hypersonic régime. 
we have essentially neglected terms of 
indeed, 


To be specific, 
0(1/.14;*) in these equations. 
sistent with the simplification introduced in Eq. (9). 
that our answer to the first 


This 1s, con- 


Hence, we may conclude 
above is affirmative. 
Probstein,'* however, the 


question In an interesting paper 


by Lees and opinion was 
expressed that Eqs. (1)—(3) are probably not valid in the 
region close to the leading edge where the shock wave 
coincides with the outer edge of the viscous layer. 
While we must admit that the present theory breaks 


down entirely in the leading-edge region, our results 


seem to indicate that the postulated type of hyper- 
sonic viscous régime may exist, at very high Mach 


numbers and a fixed Reynolds Number, over a large 


portion of the plate surface outside the leading-edge 


region (Figs. 2 and 3) where Eqs. (1)-(3) are expected 
to be valid. 

The second question above is concerning the unique- 
ness of the solutions. We know that the values of 4, 
dé/dx, and p2/p, obtained by the present method vary 
over a range with varying assumptions about the form 
of u/u%. Computations using a quadratic velocity pro- 
file,t for example, yield a value of Ay = 0.606, which is 
about two-thirds of the value used in the present paper. 
However, we found that the functional forms of the 
solutions are not affected by this choice of the quadratic 
(u/u2)-function. Furthermore, we shall point out that, 
although the present analysis is formulated with the 
explicit condition Ag # 0, nevertheless, some of our 
results have the same analytical forms as previously 
obtained by Shen.’ Therefore, it seems apprepriate to 
say that the analytical nature of the solutions is little 
affected, while their numerical values are really affected 
by the choice of the velocity distribution function. 
The quartic velocity profile has been chosen for the 
present calculations because it appears, from past ex 
perience with low-speed boundary-layer flow, that no 


1, f(O 0 and f"(0) = 


+ The boundary conditions f(1 


—A are satisfied 
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TABLE 2 a ee 
Values of Cr at Various M, for Fixed R, 
= R; = 108 ~ | - R; = 107 
M, Cr Cr Cr | M, Cr Cr Cr Cr Cr 
|Eq. (46b)] [Eq. (46a)] von Karman | [Eq. (46b)] [Eq. (46a)] [Eq. (51b)] [Eq. (5la)] von Karman 
6 0.00149 ; 0.001249 | 
7 0.001569 0.001577 0.001205 
8 0.001679 0.001690 0.001168 8 0.000420 0.000369 
9 0.001825 0.001841 0.001138 | 9 0.000428 0.000359 
10 0.002009 0.002028 0.001109 | 10 0.000441 0.000443 0.000359 
| 12 0.000479 0.000482 0. 000336 
| 16 0.000611 0.000617 0.000815 
| 20 0.000822 0.000831 0. 000299 
} 22.132 0.000966 0.000957 0.000292 
| 26 0.001164 0.001144 0.000281 
| 28 0.001268 0.001261 0.000276 
| 30 0.001372 0.001366 0. 000272 
| 32 0.001478 0.001473 0.000268 


great improvement in accuracy was gained by consid- 
ering the velocity profile as a polynomial of higher 
degree. 

It seems, therefore, that the use of the momentum 
integral method in the manner outlined in the present 
report, which has the great merit of being compar- 
atively easy, may be expected to give a fairly satisfac- 
tory estimate of the interaction effects between the 
shock wave and the boundary layer at very high Mach 
Numbers. The Mach Number range in which the pres- 
ent analysis appears reasonable must be such that 


(i) 
Xo il (C4) 7 1 
R M2 


and (11) 
M,/R,' < 1/100 


The first criterion implies that we are only considering 
the effects of such hypersonic shock waves as have a 
persistence over a length at least of the order of 1 //,°.7 
The second criterion is necessary to avoid the use of 
new boundary conditions pertaining to slip flow. 

The above discussions describe the limitations of the 
present theory. Within the hypersonic Mach Number 
range defined above, we shall calculate the variation of 
Cy with M], for fixed R, by the formulas in Eqs. (46a), 


+ The e-neighborhood of the leading edge mentioned on page 


346 has thus been given a specific definition. 
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over the entire plate surface. 














SHOCK-WAVE EFFECTS ON 
(46b), (51a), and (51b). The calculations are per- 
formed for the cases R; = 10° and 10’, with y = 1.4, 
w = 0.768, Ao = 0.91. Figs. 5 and 6 show the results 
thus obtained, which are also tabulated in Table 2. In 
these figures, Cr data on the basis of von Karman’s ap 
proximate formula’* are also given. It is clearly seen 
that, in the hypersonic speed range, the effects of the 
shock wave on the Cy of the plate certainly cannot be 
neglected. Attention should be drawn to the strong 
indication that, contrary to the conventional theory, 
which predicts a steady decrease in Cy as MM, increases, 
the present results predict a steady increase of Cy with 
M,. 

It is tempting to draw conclusions about the percent- 
age rise of skin-friction drag on a flat plate due to the 
hypersonic shock effects. However, besides the pos- 
sible inherent inaccuracy due to the assumed velocity 
profile, our numerical answers almost certainly suffer 
more or less from the various simplifying assumptions 
such as (1) Pr = 1, (2) y = 1.4, and (3) w = 0.768. 
In the large temperature range anticipated under free- 
flight conditions, none of these assumptions is expected 
to be true, and further complications arising from ef- 
fects of dissociation and ionization may also begin to 
occur. But, in the test section of a hypersonic wind 
tunnel where the perfect gas conditions are not greatly 
violated, these assumptions are satisfactory. There- 
fore, the numerical data in Figs. 5 and 6 can be regarded 
as preliminary estimates of the hypersonic shock-wave 
effects. More realistic conclusions of a quantitative 
nature probably should be deferred until more experi- 
mental facts are gathered in the Mach Number range of 
8 to 15 and higher. 

Lees and Probstein'* recently predicted a similar in- 
crease in the skin friction at hypersonic speeds. Their 
method of approach and basic assumptions are quite 
different. The 
theories in this respect is rather encouraging. 


qualitative agreement of the two 

For x» L < 1, the present theory must admit the 
coexistence of a hypersonic shock-flow régime with an 
ordinary boundary-layer flow régime over the flat plate 
(Fig. 4). The joining together of these two régimes 
must be accomplished through a transitional region. 
The present method is not suitable for exploring this 
transitional phenomena. Indeed, for the Mach Num- 
ber range corresponding to x»o/L < 1, the method of 
Lees and Probstein probably yields more dependable 
results. 

For x) L > 1, the hypersonic shock régime takes 
The interaction phe- 
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nomena between the shock wave and the boundary 
layer are believed to be more realistically represented 
by the present approach. Therefore, for the Mach 
Number range corresponding to x9/L > 1, the simpler 
method as presented in the present paper should yield 
satisfactory results. 

Examination of the shock-wave effects on the skin 
friction of a noninsulated flat plate at hypersonic speeds 
is reserved for a later paper. The corresponding effects 
on the rate of heat transfer across the plate surface 
have been shown to be rather small.!4 


REFERENCES 


‘Li, T. Y., and Nagamatsu, H. T., Shock Wave Effects on the 
Laminar Skin Friction of an Insulated Flat Plate at Hypersonic 
Speeds, GALCIT Hypersonic Wind Tunnel Report, Memoran 
dum No. 9, 1952 

* Kuerti, G., The Laminar Boundary Layer in Compressible 
Flow, Advances in Applied Mechanics II, Academic Press Inc., 
1951. 

’ Lagerstrom, P. A., Cole, J. D., and Trilling, L., Problems in 
the Theory of Viscous Compressible Fluids, GALCIT Report, 
1949. 

‘Lin, C. C., Reissner, E., and Tsien, H. S., On Two-Dimen 
sional Non-Steady Motion of a Slender Body in a Compressible 
Fluid, Jour. Math. Phys., Vol. 27, No. 3, 1948. 

5 Shen, S. F., An Estimate of Viscosity Influence on Hyper 
sonic Flow over an Insulated Wedge, Jour. Math. Phys., Vol. 31, 
No. 3, 1952. 

6 Puckett, A. E., and Stewart, H. J., The Thickness of a Shock 
Wave in Air, Quar. Appl. Math., Vol. VII, No. 4, 1950. 

7 Bethe, H. A., and Teller, E., Deviations from Thermal Equilib 
rium in Shock Waves, Ballistic Research Laboratory Report 
X-117, Aberdeen Proving Ground, 1945. 

8 Shen, S. F., On the Boundary Layer Equations in Hypersonic 
Flow, Readers’ Forum, Journal of the Aeronautical Sciences, 
Vol. 19, No. 7, p. 500, July, 1952. 

® Liepmann, H. W., The Interaction Between Boundary Layer 
and Shock Waves in Transonic Flow, Journal of the Aeronautical 
Sciences, Vol. 13, No. 12, p. 623, December, 1946. 

10 Pohlhausen, E., Zur ndherungsweisen Integration der Differ 
entialgleichungen der laminaren Grensschicht, Z. angew. Math 
Mech., Vol. 1, 1921. 

1! Goldstein, S., Modern Developments in Fluid Dynamics, 1; 
Oxford University Press, 1938. 

2 Tsien, H. S., Superaerodynamics, Mechanics of Rarefied 
Gases, Journal of the Aeronautical Sciences, Vol. 13, No. 12, 
p. 653, December, 1946. 

18 yon Karman, Th., The Problem of Resistance in Compressible 
Fluids, Fifth Convention of the Alexander Volta Foundation, 
1935. 

14 Lees, L., and Probstein, R. F., Hypersonic Viscous Flow over 
a Flat Plate, Princeton University Aeronautical Engineering Lab 
oratory Report 195, 1952. 

16 Lees, L., On the Boundary-Layer Equations in Hypersonic 
Flow and Their Approximate Solutions, Readers’ Forum, Journal 
of the Aeronautical Sciences, Vol. 20, No. 2, p. 143, February, 
1953. 








Readers’! Gorm 





RIEF REPORTS of investigations in the aeronautical sciences and discussions of papers pub- 


lished in the JOURNAL are presented in this special department. 
completed as soon as possible after receipt of the material. 
responsible for the opinions expressed by the correspondents. 


The publication will be 
The Editorial Committee does not hold itself 
Contributions should not exceed 800 


words in length. 


Further Note on the ‘‘Dynamic Stability of a 
Missile in Rolling Flight’’ 


R. E. Bolz 

Professor of Aeronautical Engineering, Case Institute of Technology, 
Cleveland 

February 24, 1953 


T° BRYSON’S COMMENTS (JOURNAL OF THE AERONAUTICAL 
ScrENCES, Readers’ Forum, January, 1953) on my paper 
“Dynamic Stability of a Missile in Rolling Flight’? (JouRNAL 
OF THE AERONAUTICAL SCIENCES, June, 1952), we believe there 
are several errors, misinterpretations, and comments that should 
be corrected. I shall discuss each of Mr. Bryson’s comments 
in the order in which he presented them. 

(1) Mr. Bryson’s difficulty in comment | can be cleared by cor- 
recting a typographical misprint in the subject paper. The 
sentence above Eq. (4) contains an incorrect word and should 
read: ‘In these equations, P; is actually the sum of the roll rate of 
the missile P about its own axis and the component of angular 
velocity w, about the Y axis due to the precessional motion re- 
sulting.’’ The meaning of Eq. (4) is easily observed by referring 
to Fig. 1 and Eq. (3) in the recent paper by Nicolaides.!. Actu- 
ally, the assumption made in neglecting w; is only that of second- 
order terms in Eqs. (3b), (3c), (8e), and (3f) and was done merely 
to identify the important Y component of angular velocity with 
the velocity normally called spin or roll in the usual gyroscopic 
axis system. 

(2) In comment 2, Mr. Bryson has unfortunately missed an 
important point. Obviously, neglecting the nonhomogeneous 
terms cannot affect the classical linearized stability analysis; this 
is well known. The inclusion of the forces and moments due to 
misaligned surfaces and jet misalignment introduces into the 
analysis the possibility of an important type of instability arising 
from the terms in the particular solution (stability criterion III) 
This is termed ‘‘resonance instability” in the more recent analysis 
of reference 1. Mr. Bryson’s inclusion of the ‘‘moment due to 
canted surfaces” in his statement is obviously incorrect 

(3) In comment 3, on the basis of consistency alone, Mr. Bry- 
son has a point. However, in actuality since nothing is known 
about the magnus forces and moments due to pitching and yaw- 
ing angular velocities and nothing can be said about these terms 
either experimentally or analytically, we felt justified in omitting 
them in the analysis in the interest of simplicity over consistency 

(4) The comment is correct. Am g is negative as defined, and, 
therefore, the sentence below Eq. (22) should have a minus sign 
preceding the term Kmg—i.e. (— King). 

(5) The first part of Mr. Bryson’s comment corrects a typo- 
graphical omission. A number 2 is missing in the denominator 
of the term in Table I, item 7, third column, Kp = Cp (.S/2d?). 
The second part of the comment is incorrect. Eq. (22) is correct 
as it appears in the original paper and needs no change. 
Kz is a normal force coefficient, it could be written in terms of the 


Since 
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usual lift curve slope and drag coefficient as 


V2 dC; pV? 
Z=—-Cp . Sa — < : 
2 da 2S 


Kia - Kia ‘+ Kp 


where Arq is proportional to the lift curve slope. If this were 
chosen to be inserted into Eq. (22), the equation could involve 


the term 
(Kre — Kp) 


rather than the incorrect one that he suggests. 

(6) Mr. Bryson’s sixth comment is also incorrect. 
stability criterion arises from the particular solution and specifi 
cally from the nonhomogeneous terms included (discussed in 
The actual importance of this third stability 


The third 


second comment). 
criterion we had failed to notice in the analysis and is clearly em 
phasized in the recent reference 1. 

(7) In comment 7, Mr. Bryson again failed to note an im- 
portant point. It is not the transformation from time to dis 
tance along the flight path that has anything whatsoever to do 
with the conclusion under #7 on page 401 of the reference paper 
namely, that dynamic stability is independent of velocity (except 
for variation of the coefficients with Mach Number). The point 
that Mr. Bryson apparently missed is that the analysis is valid for 
variable linear velocity V and variable roll velocity P = P(V) 
(but small roll acceleration). The term P; = P/V in the equa 
tion is then, in general, a variable. However, according to the 
discussion on the bottom of page 398, the point is emphasized 
that, over the Mach Number range from 1.3 to 4.0, P, has been 
shown to be essentially constant? thus rendering the equations as 
written independent of the velocity V and justifying my conclu 
sion. The author has no data concerning the effect of velocity on 
P, in the high subsonic Mach Number range, thus prohibiting the 
conclusion from being stated for this range 

(8) Mr. Bryson’s eighth comment is one of personal preference 
Using the usual reference axes system that he mentions will elim- 
inate several terms from Eqs. (8a) and (8b). Subsequent equa 
tions, however, are of identical complexity in both cases, and the 
question becomes one only of preference. 

(9) Mr. Bryson’s comment of item 3, Tables I and II, is grate- 
fully acknowledged. The statements were intended and should 
“Forces (moments) due to rolling combined with angles 
The second part of Mr. Bryson's comment 


read: 
of pitch and yaw.” 
Items (3) and (4) in Tables I and II are not at all 
the same or equal. A,, and the moment A», arise from wing- 
and fin-tip suction forces resulting from the rolling velocity. As 
such, they are additive to the magnus coefficients arising from 
the rolling fuselage. The K,g and Ama terms are damping de- 
rivatives arising from the surface potential on the lifting surfaces 


is incorrect. 


for accelerated motion and are essentially independent of the 


rolling motion of the missile. 
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(10) I concur with the part of Mr. Bryson’s comment con- 
cerning a vagueness in conclusion (D) on page 400 of the original 
paper, except for his underestimation of the size and influence 
of magnus moment. For spin stabilized missiles, this factor ap- 
pears to be of prime importance in its contributions to dynamic 
instability. Actually, referring back to Eq. (24) of the original 
paper, conclusion (D) should be clarified as follows: 

(D) In general, a statically unstable rolling missile can be made 
dynamically stable by high roll rates, although this method of ob- 
taining stability is of practical interest primarily for the wingless 
and tailless spinning body. For a statically stable missile, dy 
namic instability can be induced by high roll rates because of 
the important influence of the magnus moment, which from 
experimental indications is not at all negligible in many cases. 

(11) It should be added that Mr. Bryson’s analysis, first, is not 
free from error and, second, misses the elements intended in the 
original paper—namely, the influence of flight and rolling velocity 
variations (except on the coefficients themselves), thrust, canted 
fins, and assymmetry conditions.! The errors in Mr. Bryson’s 
analysis are mathematical errors in the equations for b, and the 
stability condition S.. The magnus force should not appear in 
any of these equations. 


REFERENCES 
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Turbulent Channel Flow from a Variational 
Principle 


Stanley Corrsin 
Associate Professor of Aeronautics, Department of Aeronautic 
The Johns Hopkins University, Baltimore 


February 2, 1953 


= B. MILLIKAN’S DEMONSTRATION that the steady parallel 
. laminar flow solutions of the Navier-Stokes equations are 
those giving minimum total dissipation! can be easily modified 
to include corresponding steady-state turbulent flows. 

Consider, for example, fully developed plane flow between 


parallel walls. The variational condition is 


Ph/2 
5 ody =0 (1) 
h/2 E 
where 
= r(dU/dy) (2) 


is the total rate of dissipation in (laminar or turbulent) steady 
state flow and 7 is the total mean shear stress. U(y) is mean 
velocity and y is spatial coordinate perpendicular to U 

The “constraint” is conservation of mass flow rate, 

 _— 
f . U dy = constant (3) 


9 
. < 


Carry out the variation with respect U(y) holding r(y) fixed.* 
The Euler differential equation corresponding to Eqs. (1) and (3) 


1S 


A’ — (dr/dy) = 0 (4) 


where \ is the Lagrange multiplier 


* Holding 7(y) fixed in the variation process is the modification referred to 
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Putting \ = OP/Odx, the mean longitudinal pressure gradient, 
we obtain the differential equation for either laminar or turbulent 


flow between parallel walls. 
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The Gas Bending Stress on a Blade in an Axial- 
Flow Compressor 


John Mitchell Stephenson 

Research Associate, Division of Engineering, Brown University, 
Providence, R.I. 

February 16, 1953 


(1) SuMMARY AND CONCLUSIONS 


A“ APPROXIMATE EXPRESSION is found for the maximum bend- 
ing stress imposed on a compressor blade in terms of the 
forces applied by the gas stream and thence in terms of the mass 
flow, work done, and blade geometry. For a given aerodynamic 
design, which implies a given solidity, the following results ap 
pear: . 

(a) The maximum bending stress depends on the thickness of 
the blade root section but is independent of the chord. 

(b) The weight of the blade row is directly proportional to the 
blade thickness and is also independent of the chord. 

(c) If the blades are hollow, the weight is reduced, but the 
stress may be increased. However, a figure of merit proportional 
to the ratio bending force/blade weight (for a given maximum 
stress) is improved by a factor of the order three, irrespective of 
the wall thickness, provided this is small. 

(d) The natural bending frequency of a blade is independent of 
the chord and depends only slightly on the material. 

Three conclusions can be drawn: 

(a) For a given stage aerodynamic design, a reduction of blade 
weight can only be achieved by raising the allowable gas bending 
stress level or by using hollow blades. (However, the selection 
of the aerodynamic design parameters for minimum weight, as- 
suming a given stress level, is a much more difficult problem. ) 

(b) Since the weight of the blades is a large fraction of the total 
weight of a stage, especially for small hub-tip ratios, a designer 
only incurs a small weight penalty by using extremely large blade 
chords, while two advantages accrue. First, he increases the 
Reynolds Number. Next, by reducing the thickness/chord 
ratio in this way he is able to use higher gas Mach Numbers 
The chief disadvantage is the increased axial length of the shaft. 

(c) The only effective ways of raising the natural bending fre 
quency of a blade are to thicken it or to make it hollow, which in 
creases the frequency by a factor of V 3, again irrespective of the 
wall thickness 

It is frequently supposed that large chords 
values of the blade height/chord ratio—imply low stage effi 
ciencies by analogy with aircraft wings 
given that the drag of the secondary vortices depends on the 
chord, and the total blade area exposed to gas friction depends 
only on the blade height, hub/tip ratio, and solidity. The loss of 
efficiency due to tip clearance is also independent of the chord 


and, hence, low 


But no proof has been 


Only the wall drag certainly increases, and this accounts for less 
than a third of the total loss. 

The presence of centrifugal loading does not alter the analysis 
given here but merely lowers the level of the maximum allowable 


bending stress for any given material 
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Fic. 1. 


System of coordinates. 


(2) ANALYSIS: THE STRESS 


To discuss the moments acting on a blade fixed at one end, 
we use the (left-hand) system of coordinates shown in Fig. 1, 
taking the origin at the centroid of the root section. The cen 
troids of the other blade sections do not then necessarily lie on a 
radial line or even on a straight line, but the effect of this will be 
neglected. 

The aerodynamic force components on any section of the blade 
are denoted as in Fig. 2. Since the principal axes of the blade 
sections are always nearly parallel and normal to the chord lines, 
we calculate the bending moments about the axes 2’ and y’, in- 
clined at the root section stagger angle ¢ to the axial axis z and the 
tangential axis y; thus, 


My = fo" Y'-x dx 
My = ,* Z'-x dx 
The stress at any point in the root section is, therefore, 
o = My-2'/Iy — Mz-y'/Ts 


where the /’s are the second moments of area and the positive 
Now Jy’ is of the order 100/,’, 
and the force Z’ is usually small compared with Y’. 


sign in o indicates tensile stress. 
Therefore we 
may write, approximately, 


o0 = —Mz/uc3 


where wc* = J,’/y’ is the modulus. We are only concerned here 

with the value of uw, corresponding to maximum stress, which js 

usually found near the leading and trailing edges. Thus, 
; heap , ee 

—opc’ = So "(Y cos ¢ + Z sin ¢)x dx 


If we write Z/Y as 8 and assume that, approximately, 8 is inde. 
pendent of x, we have 


—ouc*® = (cos ¢ + B sin ¢)- , Vx dx 
Now if the average work done is CpAT per unit mass flow, we have 
Y-Buh = Q-CpAT 


where B is the number of blades, u is the mean blade speed of 
rotation, and Q is the mass flow. Therefore the stress equation 


becomes 
—opc? = (cos ¢ + Bsin £)-(Q-CpAT/2u)-(h/ Be) 


Note that Bc depends on the solidity of the blades, which we may 
assume, like all the other terms on the right side, to be fixed by 


the aerodynamic design of astage. Therefore we may write 


where F is a force proportional to the bending force on the blade, 
independent of the size or shape of the blade root profile 


(3) THE Mopvutus 


In general, uw can only be found by analysis of the particular 
profile shape used for the blade root section; but toa first approxi- 
mation we may write, fora solid blade, J,’ « ct? and y’ « t. Hence, 
uc'« cf?, and we find that up & k(t/c)?, where & is a number of the 
order ().1 and ¢ is the maximum thickness of the section. (This 
formula is exactly true for an ellipse, where k = 7/32, and fora 
rectangle, where k = 1/6). 

Thus we have, approximately, 


This is our first result-—that the stress is independent of the chord 
For a hollow blade, made of sheet of a constant thickness d, 


we have approximately, if d< ¢, /,,« ct?d and y’ = ?; 


and we 











Fic. 2. The force component diagram 
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find that uc? ~ 6ktd, where k is the same number as for the solid 
blade. (Again, this formula is exact for an ellipse and for a rec- 
tangle, if d<t). Thus the stress formula for hollow blades is 


—6otd = F/k 


Hence, if the wall thickness is less than ¢/6, the maximum stress 


is increased by hollowing a blade. 


(4) AN EXAMPLE 


The application of the approximate stress equation may per- 
haps be made clearer by a numerical example. Consider a com- 
pressor Stage with the following properties: Q = 200 lbs. per sec. ; 
Cp = 10,900 ft.?/sec.; AT = 20°K.; u = 750 ft. per sec.; h = 10 
in, Bc = 80 in.; and assume that k = 0.1 and 

(cos ¢ + Bsin ¢) = 1 
Then we have F = 113 Ibs., so that 
—ol? = F/k = 1,130 Ibs. 
Therefore we may draw up a table (Table 1) to show the greatest 
number of solid blades which can be used for a given value of the 


root thickness/chord ratio and of the maximum allowable 


bending stress. Therefore, using a smaller number of blades 


(of larger chord ) reduces either ¢/c or the maximum bending stress. 


TABLE | 


Max. o 
1,000 Min. ¢, t/c at Root 
Lbs. per Sq.In In. 0.08 0.10 0.12 
20 0.238 27 33 40) 
24 0.217 29 37 44 
28 0.201 32 40 48 


(5) THe WEIGHT 
If blades always have the same height and taper ratios of thick 
ness and chord, we can say that the weight W of the blades in a 


row is proportional to the blade root area A: 


We« BA 
Now for a solid blade, A ™ act and, for a hollow blade, A ™ 2acd, 
where a is a number of the order 3/4. (For an ellipse, a = 7/4; 


fora rectangle,a = 1). Therefore, since, again, c only occurs in 
the combination Bc, we have our second result that the blade 
weight is independent of the chord. 

Now we define a dimensionless figure of merit F/W, which is 
proportional to the ratio of blade loading /weight; we require this 
to be as high as possible for any given value of the stress ¢. But 
we have seen that 

F = — ou 2 
Thus F/W« pc?/BA = A,say. Using the values of u from para- 
graph (3), we have 

A = kt/aBc for the solid blade 

A = 3kt/aBc for the hollow blade 


| 


Therefore, with this criterion, a hollow blade is about three times 
better than one externally similar but solid. 


(6) VIBRATION 


The alternating gas forces on a blade may be assumed to pro- 
duce stresses that are approximately proportional to the steady 
gas bending stresses, so that they are also independent of the 
chord. Moreover, the lowest natural frequency of a blade in 
bending is f, where 


f? ~ El,z’/yAh* 


where E is Young’s modulus and y is the density of the blade 
material Hence, 


f? ~ Epc3t/y Ah* = EABct/yh' 
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or 
f?« EAt/y 


which is independent of the chord (our third result ) 

Thus, there are three ways of raising the natural bending fre 
quency (to avoid resonance, for example): 

(a) Increasing the thickness ¢. 

(b) Using hollow blades, with their large value of A {see para 
graph (5)]. 

(c) Raising the value of E/y by changing the blade material 
However, E/y is roughly the same for many metals, and this is 
rarely a useful move 

I should like to thank Prof. Max Schilhansl for his interest and 
for drawing my attention to his paper on this subject.! 
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Additional Comments on ‘“‘A Study of 
Transonic Gas Dynamics by the Hydraulic 
Analogy’”’ 


J. R. Bruman 


Jet Propulsion Laboratory, 
Pasadena, Calif 


January 27, 1953 


California Institute of Technology, 


— PART OF THE SUBJECT PAPER! describing results obtained 
by the writer? and the subsequent discussion in the Readers’ 
Forum* * appear to foster a slight misunderstanding based on 
ambiguous use of the term ‘‘wave.”’ 

The free-surface wave under discussion is a continuous, peri 
odic, frictionless disturbance not related to the shock wave or 
its free-surface analog, the hydraulic jump. The only connection 
between the two phenomena is the fact that the minimum prop- 
agation velocity for the jump equals the value approached by the 
surface wave as its wave length increases. 

Since experiments with the hydraulic analogy are usually de 
signed to represent shock waves by hydraulic jumps, it is desir 
able to avoid disturbances propagated at velocities other than 
that of the jump being observed. It is immediately seen that 
spurious disturbances of sufficient wave length will at least be 
propagated at the velocity of the slowest hydraulic jump. The 
heuristic notion that this might be the most desirable test condi 
tion led the writer to propose and investigate the water depth 
for which the wave length required is a minimum.” For water, 
this depth is about 0.225 in. in theory—and proved, in fact, to be 


for observation of hydraulic jumps. However, it 


the optimum 
seems difficult to reach this conclusion about shock waves, repre 
sented by hydraulic jumps, solely from a consideration of periodic 
surface waves and their velocities 

In the case of those observers® ® who chose to apply the anal 
ogy with a moving stream and stationary model, the !/,-in 
depth was impractical, but this did not prevent successful work 


at larger depths with sufficiently large models 
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Location of Detached Shock Waves Ahead of 
Plane Bodies 


W. E. Moeckel 
Lewis Flight Propulsion Laboratory, N.A.C.A., Cleveland 
January 29, 1953 


ia A RECENT READERS’ FORUM LETTER, Professor Truitt! pro 
posed what he believed to be a simple new method for pre- 
dicting the location of detached shock waves ahead of wedges in 
transonic flow. Examination of his equations shows, however, 
that his method is identical to one of the methods (the ‘‘geomet- 
rical method’’) proposed in reference 2._ Thus, in reference 2, the 
equation for shock location is given as (in the notation of refer 
ence 1) 

(S + L)/t = cot 0’ (1) 


where (.S +- L) is the distance from the shock to the body sonic 
point, / is the ordinate of the sonic point, and 6’ is the maximum 
wedge half-angle for an attached shock at the prescribed Mach 
Number. Eq. (1) is exactly the equation from which Professor 
Truitt derives his final expression. 

Professor Truitt appears to have misinterpreted an assumption 
made in reference 2 and concludes that the misinterpreted as- 
sumption is invalid. The assumption made was that the axial 
distance from the body sonic point to the detached shock wave is 
independent of nose shape. This assumption is embodied in 
Eq. (1). Since the methods of reference 2 are based on this as- 
sumption and since Professor Truitt’s method is identical to one 
of those methods, the good agreement between experiment and 
theory shown in reference 1 confirms, rather than disproves, the 
assumption made in reference 2. Professor Truitt’s conclusion 
that there is an effect of nose shape on shock location results from 
his use of distance from wedge vertex instead of distance from 
wedge sonic point in plotting the data. 
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Comments on a Method for Location of 
Detached Shock Waves Ahead of Plane Bodies 


E. V. Laitone 
Associate Professor, University of California at Berkeley 
January 29, 1953 


ia REFERENCE 1, TRUITT ERRONEOUSLY BELIEVES that he has 
found a new equation for predicting the location of a detached 


shock, one that includes the effect of nose shape. However, all 


he has done is to replot Moeckel’s geometric result of reference 


AERONAUTICAL 


MAY, 


1953 


SCIENCES 


2 toa different scale. This may be shown by writing Moeckel’s 


Eq. (2) of reference 2 as 


, t t 
tan?’ = : - tan@ = | 
( + 2 (7) 


corresponding to Fig. 1 of reference 1, where @ is the wedge half 
angle and @’ is the corresponding maximum angle for which the 
shock becomes detached. Then, solving for S, the detachmen; 
distance from the nose of a wedge having total length L (see Fig. |, 
reference 1), we obtain 

S/L = (tan 6/tan 0’) — 1 (2 


which is identically Truitt’s Eq. (2) in reference 1. This does 
not actually predict any direct effect of nose shape, as may be 
seen by the original form of Eq. (1), using the detachment dis 
tance measured from the sonic line (S + ZL). These same equa 
tions are also applied by Moeckel to bodies of revolution, simply 
defining 6’ as the cone half-angle for which the shock becomes de 
tached. It is obvious that Eq. (1) is much more general thar 
Eq. (2), which can only be applied to a wedge or cone 

Actually, these results are only valid for relatively blunt bodies 
at the higher transonic Mach Numbers as shown in reference 3 
The only results valid for very slender bodies at Mach Numbers 
near unity are those obtained by the small perturbation method, 
which I originally derived in 1945 (see reference 4) and later used 
in references 3 and 5. This small perturbation method gives the 
following equation for the conditions shown in Fig. 1 of refer 
ence |: 


Z 4\/2 47 
Infl+—j]= ( 
S (vy + 1)0 
As shown in reference 3, Eq. (3) does predict the correct effect 
of nose angle but only for very slender wedges (@ < 5°). Also 
note that Eq. (3) satisfies the transonic similarity relations 


Wo — 1) 3 
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Note on Effects of Density Fluctuations on the 
Turbulent Skin Friction of an Insulated Flat 
Plate at High Supersonic Speeds 


E. R. Van Driest 

Aerophysics Department, North American Aviation, Inc., 
Calif 

February 2, 1953 


Downey, 


= THEIR READERS’ FORUM REPORT, “Effects of Density Fluctu- 
ation on the Turbulent Skin Friction of an Insulated Flat Plate 
at High Supersonic Speeds” (JOURNAL OF THE AERONAUTICAL 
SCIENCES, October, 1951), T. Y. Li and H. T 
that the density fluctuations occurring in a compressible turbulent 


Nagamatsu state 


boundary layer at high speeds are not adequately taken into 
account in the analyses presented by certain authors, including 
that of the writer.! The purpose of the present note is to dem- 


onstrate that the statement of Li and Nagamatsu is not justi 
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READERS’ 


fied, because the procedure followed by them appears to be phys- 
ically inconsistent and, therefore, may lead to results that can- 
t be expected to agree with experiment. 


no 
In the writer’s paper,' the following equations for continuity, 


momentum, and energy were derived to describe the mean 
steady flow in a thin turbulent boundary layer: 
(0/Ox) (pu) + (0/Oy) (py) = O (1) 
Ou Ou i - 
ou + pv = [— (pv)’u’] (2) 
Ox oy oy 
oT oi re) 7 — On t 
pu T pv = — (pv)? — (pv) u (3) 
Ox oy oy oy 


where p is the density, u is the velocity component in the direc 
tion x parallel to the plate, v is the velocity component in the 
direction y perpendicular to the plate, and 7 is the enthalpy. 
The bars indicate temporal mean values and the primes the in- 
stantaneous fluctuation from the mean. It is seen that the mo 
mentum and enthalpy equations, Eqs. (2) and (3), are consistent 


with the continuity equation, Eq. (1)—i.e., the equations are 


consistent as far as the averaging process is concerned. The 
terms 
mae 
(pv) uo =r (4) 
rr = 
—(pvj)jt =@ (90) 


thus may be properly called the Reynolds stress and Reynolds 


heat transfer, respectively. If, now, one puts 


(pv)'u’ = (Ou /dy) 
(pv)’t’ = w(Ot/dy) 
and assumes that 07/O0x = Oand ¢ = « (turbulence Prandtl Num 
ber = 1), it will immediately follow that 
t+ (u?/2) = constant (6) 


for an insulated plate 
It is noted that density fluctuations are still included in Eqs 
1-6) and, in particular, in expressions (4) and (5) 


Expression (4) is next expanded; whence, 


r = —(pv)'u’ = — pu'v'’— 3d p'u' p'u'y’ (7) 
which reduces to 
rT = —pu’d’ p'u'v’ (8) 
when * = 0 fora thin boundary layer, and to 
T= pu ’y’ (9) 


when the triple correlation p’u’'v’ is neglected. Thus the density 
fluctuation appears only in the triple correlation term p’u‘v’, 
and, therefore, if p’u’v’ is small compared to p u’v’, Eq. (9) ade 
quately describes the compressible flow Reynolds stress. Since 
Eq. (9) was used to derive a friction law in reference 1, then it is 
true that density fluctuations were neglected by the writer in his 
final analysis; however, not to the extent indicated by Li and 
Nagamatsu 

Eq. (1) of the note by Li and Nagamatsu, on the other hand, 
results from adding to the above Eq. (7) various degrees of the 
term @ p’v’, which does not even appear in the writer’s analysis 
In their work? the term @ pv’ came about as a result of using 
p-% and p-d for mass flux, rather than pu and py as required by 


continuity [see Eq. (1) above]. If a p’v’ is added to Eq. (7), 
then 
r= pu'v’ ui p'v’ D pu’ p'u'v’ (10) 
Hence, when 5 = 0, 
r= pu'v’ ui p’v’ p'u'y’ (11) 


Which is equivalent to 
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in contrast to the writer’s —(pv)’u’ 

The difference in physical interpretation of (pu )'v’ and ( pv’ )u’ 
appears clear, the latter being consistent with the transverse 
momentum-transfer concept. It is because of this consistency 
with regard to the transport (flux) concept that the writer be 
lieves that Eqs. (1)-(3) are properly formulated 

If, now, the mixture concept be applied to both momentum and 


velocity gradients and if, further, p’u’ be neglected, then 


r = —/? [d(p-u)/dy|-(di/dy) (13) 
which is the same as the Li- Nagamatsu expression for a = 1 in 
their note. When it is assumed that / = ky and a mean friction 
law is derived following the procedure in reference 1, there re 
sults, for an insulated plate, * 


0.242 (1 — A?)/? @(A) ’ 
; . = logy (Re-Cy) + 


Cy d 
0.76 logio (1 A?) + = logw (1 + A?) (14) 
where 
(A 2 tanh”! sl tanh”! 2 

QA) = VS ta at 
Vi+. Vi+» 

l(y — 1)/2|M .? 

? = 
1 + [(y 1)/2|.M 


In Fig. 1 is plotted the ratio of the compressible-flow skin 
friction coefficient to the incompressible-flow skin-friction coeffi 
cient, Cy,/Cy;, obtained from Eq. (14). Also plotted are the anal 
yses of Wilson‘ and the writer,' both based on Eq. (9 A strong 
divergence of Eq. (14) from the two theories using Eq. (9) is evi 
dent. (The difference between Wilson’s analysis and that of the 
writer is that Wilson‘ used the similarity law for mixing length, 
whereas the writer used the linear relation / = ky.) The fact 
that the data‘ fall near the theories of Wilson and the writer 
should be sufficient evidence that the use of Eq. (12) instead of 
Eq. (4) may lead toerroneous conclusions. To show this result was 


the purpose of the writer’s comments in reference 3. It never was 
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expected by the writer that Eq. (14) would be seriously considered, 
because of the physical difference between (pv)’ u’ and (pu)’v’. 

It should be pointed out further that including the term 
ii p'v’ not only brings about a definite change in the shear equa- 
tion from + = — pl*(dui/dy)* to Eq. (13) but also leads to diffi- 
culty in the derivation of the energy relation Eq. (6). In the 
opinion of the writer, either one can disregard density fluctu- 
ations altogether in accordance with Crocco,’ or one should follow 
the course suggested above by the writer. Liand Nagamatsu did 
not derive the energy relation in their note or main report ;? how- 
ever, it is the writer’s belief that, if they had attempted to derive 
the energy relation using their interpretation of density fluctu- 
ations, they would have encountered considerable difficulty. 

Summarizing, density fluctuations have been neglected in the 
derivation of Eq. (9), but only to the extent of neglecting p’u'v’ 
and not a p’v’, because a p’v’ does not even appear in the writer’s 
analysis. At high Mach Numbers, p’u’v’ may become important. 
However, there are many other approximations irvolved in the 
turbulence theory all the way down the line which may also 
become increasingly important; for example, in the derivation of 
friction laws, it is generally assumed that p ~ 1/T to the neglect 
of p’7’, where T is temperature. All in all, at least it seems 
clear that certain authors, including the writer, have adequately 
taken density fluctuations into account by use of Eq. (9) and, 
therefore, that the referred-to statement by Li and Nagamatsu is 
not justified. 
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Further Comments on ‘‘Aerodynamic 
Coefficients of an Oscillating Airfoil in Two- 
Dimensional Subsonic Flow’’ 


P. F. Jordan 
Royal Aircraft Establishment,* Farnborough, Hants, England 
February 27, 1953 


IT THE SUBJECT PAPER,! Timman, van de Vooren, and Greidanus 
published numerical values; these values differ from the 


standard values—i.e., those published by Dietze.? Fettis’ re- 
marked that the new values! must be erroneous because they do 


not, for \/ > 0, satisfy the relation f 
ka = tw(ma + kp) (1) 


which is, however, satisfied by the standard coefficients.2. In a 
recent letter,4 the authors of the subject paper! maintain that 
their values must be correct; they state among other things: 

* Acknowledgment is made to the Chief Scientist, British Ministry of 
Supply, for permission to publish this note. 

+ The notation, identical to that of reference 1, is based on the Kiissner 


notation, but the reference point is the mid-chord. Fettis’ derivation of 


Eq. (1) is not yet available to the writer. 
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(a) The relation (1) holds only for incompressible flow. 

(8) The standard results? may well be erroneous for J/ > 0 bp. 
cause Dietze used an approximate method. 

(y) Dietze’s iteration method does not converge if both M and 
w exceed 0.7. 

The writer has the following comments to make: 

(a) The relation (1), taken by itself, isnot new. For Jf > Lit 
was stated independently by Temple and Jahn’® and by the 
writer;® it was then observed® that Eq. (1) is valid for Mf = 0 byt 
has no practical significance there. 

(b) The relation (1) can be proved for any Mach Number J 50 
that reference 1 must contain an error. Let y, be the (nondimen- 
sional ) lift distribution due to the downwash w,, = x” at the fre- 
quency given; then Eq. (1) is easily shown to be equivalent to 


Py (v1 + xy) dx . (woyi + wiyo)dx =0 (2 


The proof of Eq. (2) for M > 1 is elementary; it follows from the 
form of the basic equation 
"Ez 


y=wt+ So w(x — &)R(E) dé (3 


by means of partial integration. The proof for 17 < 1 is more 
involved but can be obtained by modern reverse-flow consider- 
ations. In fact, Eq. (2) is readily seen to be a special case of the 


general theorem 
ye (wy — wy) dS = 0 (4 


see Flax’, Eq. (29). Eq. (4) is valid for any plan form S aid any 
two downwash distributions w and iw; 7 arises from @ in reverse 
flow. [Flax’ also states a specific relation equivalent to Eq 
(1).] 

(c) The argument that leads to the statement (8) in reference 4 
is hard to follow. Not only Dietze but a number of other au- 
thors, notably Schade,’ have calculated the aerodynamic coeffi- 
cients in question, using approximations of essentially different 
forms to the exact kernel of the Possio equation. However, these 
different approximations are so good that they should not lead 
to appreciable numerical errors; in fact, the numerical values ob 
tained by different authors agree well enough. 

(d) The Possio equation is, of course, equivalent{ to the 
linearized differential equation on which the theory of reference | 
is based directly. Recently, Jones? has recalculated (for w < 
0.5) the coefficients in question by means of a simple method* ! 
based on the transformed differential equation; his results con 
firm not only Dietze’s values? (17 < 0.7) but also Schade’s values* 
for M = 0.8.** (The latter values had not, to the writer's 
knowledge, been confirmed before. ) 

(e) The writer is surprised by statement (y). For (\/,w) = 
(0.7, 0.5) the ratio of consecutive corrections is still ~0.5 in 
Dietze’s iteration. Hitherto the writer had the impression that 
Dietze’s iteration would be convergent for any set (1/7 < 1, 


This impression was based on the following fact: the relation 


Yn = sin ng if w, = V1 — M* cos ng (5 


which is correct for w = 0, is for w # 0 an approximation that 
becomes increasingly close as increases for a given set (M,w) of 
parameters. 

(f) However that may be, the region in the M],w-plane in 
which Dietze’s method may be applied in practice is certainly 
Hence, and not 


only for this reason, the error that has occurred in the elegant 


smaller than it is for the method of reference 1. 


analysis of reference 1 is much to be regretted, and the usefulness 
of the “reverse-flow relation,” as Eq. (1) should properly be 
called, becomes apparent in that it should enable the authors of 
reference | to trace that error. 


t It can be verified for instance for M = 1, where the general solution of 
the differential equation exists, that this solution also satisfies the Possie 
equation 


** Jones’ has also calculated the coefficients for M 0.9, aw < 0.3; these 


coefficients again satisfy Eq. (1) 
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Remarks on Thermal Characteristics of 
Boundary Layers* 


Martin Bloom 

Assistant Professor of Aeronautical Engineering, Polytechnic 
Institute of Brooklyn 

January 30, 1952 


NOTATION 


a = b/7 
b = empirical constant (5,526°R. for air) 
cf = local friction coefficient (as in reference 2) 


heat capacities at constant pressure and constant volume 


(eT 
h enthalpy (| Cp ar) 
vo 


he local heat-transfer coefficient (as in reference 2) 
k = thermal conductivity 
\ local Nusselt Number (/-x/k) 
r Prandtl Number (yc,)/k) 
R = perfect gas constant 
R; = local Reynolds Number (pux/p) 
rT = temperature recovery factor (7, T1)/To Ti) 
Th = enthalpy recovery factor (he — hi) /(ho — hi) or 
(°T ¢ / To 
J, ae f. cp a7 
1 = absolute temperature 
u = streamwise velocity 
x streamwise coordinate 
p mass density 
“ = absolute viscosity 
Cp/cr 
Subcripts: 
e adiabatic surface 
] local values at edge of boundary layer 
0 isentropic stagnation conditions 
u value at surface 


DISCUSSION 


_ EXACT CALCULATIONS of laminar boundary-layer character 
istics appearing in references 1-3 and the recent reports of 
carefully measured temperature recovery factors (e.g., references 
4and 5) prompt the following remarks: 

(1) It bears restating that the variation of P with 7 for air 
has not been closely established because of uncertainty in the 
quantitative variation of k with 7. Data used in references 1-38, 
based on references 6-9, indicate a decrease in P from about 0.74 
at 500°R. to about 0.63 at 6,000°R. On the other hand, lower 
values of k reported in references 10-12 (shown in Fig. 19, refer- 
ence 1 and Fig. 3, reference 13) would lead to higher values of P 
and give credence to Knuth’s curve (reference 14, based on refer- 
ence 15), which shows P increasing slightly from about 0.74 at 
540°R. to about 0.78 at 5,400°R. An increasing trend is also indi- 
cated in references 16 and 17 (since y is reliably known to decrease 

* Done in connection with U.S.A.F. contract AF 33(616)-116, project 
52-670A-86, Wright Air Development Center 
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with increasing 7). The estimated decrease in P, shown in 
Fig. 5, reference 2, is based upon the assumption that k/yc, does 
not vary with 7. This differs from the Eucken assumption 
(that k/ucpy ~ ~ 9y — 5), which is used in references 14 and 15, 
and leads to an increase in P. Young and Janssen! noted the un- 
certainty involved in evaluating P precisely but concluded that 
the effect of the variation was minor in their application. Thus, 
it seems to be most reasonable, considering the present state of 
the art, to treat P as a constant, even in accurate calculations. 
This gives added emphasis to the important work of Crocco'’ who 
used the Sutherland 4 — T formula, treated cp as an arbitrary 
function of 7, accurately known, and maintained P constant 
Furthermore, attempts to correlate data on rr with functions of 
P evaluated at various temperatures, such as 7,, 7», or 7), must 
be viewed with care. This is particularly true if data are ob- 
tained from tests in which the temperatures produced are much 
lower than those encountered in practice. 

The following commonly referred to values are given for con- 


venient comparison: 


P P > P' i 

Approximate low-temperature value 0.74 0.86 0.91 
Possible high-temperature values 0.63 0.79 0.86 
0.78 0.88 0.92 


(2) One of Crocco’s significant results, derived for laminar 
flow on a flat surface [Eq. (40), reference 18] is given as follows: 


P'/? = (he — In)/(lo hy) rh (1) 
where hy is defined by the relation 
ho — hy = u,?/2 (2) 
Values of rr calculated from Eq. (1), by assuming 7, = 0.85 
and 7, = 400°R., differ from r;, by less than 2 per cent for 1; - 


10. Thus, values of 7, can be calculated closely from the ex- 
pression 
(Te — Ti) MT — Th) = rr Sr = PP (3) 
where 7) must be calculated with cp as a function of 7. In prin- 
ciple, however, rr is seen to be slightly dependent upon .\/; and 
7, for laminar flow on a flat surface, assuming P is constant. 
Since the error incurred by neglecting this variation is probably 
of the same order as the uncertainty concerning the value of 
P’’?, it is again seen that the value of precise correlations of meas- 
ured values of rr with functions of P is questionable at present 
In fact, an experimental investigation of the validity of Eq. (1), 
which would require the simulation of elevated temperatures, 
might provide information on the variation of P with 7 
From a computational standpoint, Eq. (1) can be used as 
readily as Eq. (3). Moreover, the accurate analytic expressions 
for cp, h, and Ty given by Donaldson" can be used to advantage 


in connection with Eq. (1). They are 


Cp/R = (7 2) + a%e"(e* — 1)~? (4) 
h/R = (7/2) Ti) + (2/7)a(e* — 1)7'J (5) 
To/T, = 1+ (y1 7)M,? — (2 7 Ja, [(e — 1)! 
(e™ — 1) 1] (6) 
Analogously, Eq. (1) yields 
T./T, = 1 + ri(yi/7)Mi2 — (2/7)ai[(e” 1)-' - 
(e”! I-44 () 


For values of 7) 660°R., corresponding to altitudes less than 
I g 


4X 105 ft., Eqs. (4)-(7) give 


cp,/R = 7/2, ym = 7/5, In = cp,Ti (8) 

T/T, = 1 + (M,2/5) — (2/7)ay(e* — 1)! (9) 

T/T, = 1 + 7i( M,2/5) — (2/7)ai (e** 1)-! (10) 

It can be noted that the curve of 7, 7, shown in Fig. 13 of 


reference 2, for no dissociation, was calculated by using an equa 
tion identical to Eq. 1 given here, with a somewhat arbitrarily 
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assumed constant value of r, = 0.85. It is of interest to com- 
pare the values shown in Fig. 13 of reference 2 with those that 
were calculated precisely by the numerical procedure and given 
in Table 1 and in Fig. 16A of reference 2. For values of u; X 
10~% (which is close to 7, in the case treated) less than 8, the 
values given by Table 1? agree with those shown in Fig. 13.% 
But for uw; X 10°* = 9.8 and 12.5, the values given by Table 1? 
are about 1,900° and 6,100°R. higher, respectively, than those 
plotted in Fig. 13;? while, for «4; X 10~* = 12.4, the value given 
by Fig. 16A? is about 900°R. lower than that of Fig. 13.2. The 
value of 7, — 7, determined from Table 3, reference 1, for M, = 
12.4 agrees more closely with those of Figs. 13 and 16A? than 
with that of Table 1.2. For JJ, > 8 the values of 7, given by 
Table 3! are not in close agreement with those of Table 1.2. The 
significance of these differences is unexplained. 

It should also be noted that the temperature recovery factor, as 
defined in reference 1, is based on 7) determined with c, con- 
stant. Its decrease with 1/7, would be less pronounced if it were 
defined according to Eq. (3), here. For example, for M, = 
9.59, 7, = 388°R., and 7, = 5,135°R. (Table 3'), the value rz = 
0.79 is obtained from Eq. (3), as compared with 0.66 given in 
reference 1. 

In connection with reference 2, it is interesting to observe that 
the effect of dissociation on the rate of heat transfer is small, at 
least for 4, < 13, if the surface temperature is less than about 
3,500°R. This is true because of the counterbalancing effects of 
dissociation on the heat-transfer coefficient and the temperature 
potential which, individually, are appreciably affected by disso 
ciation. The effectiveness of the concept of a heat-transfer co 
efficient, therefore, appears to be reduced under these circum- 
stances, as it is when the surface temperature is variable. 

(3) Crocco also gave the following expression for the heat 
transfer coefficient [Eq. (41), reference 18]: 

Tors? 
(11) 


(NP—'/R,—'/2),/(0.5¢;Rr'/*)1 = (lw — hea" Te 


He noted that Eq. (11) is virtually uninfluenced by the viscosity 
law, despite the fact that cy depends strongly on the viscosity 
law. However, cy is not greatly affected by cp and P. There- 
fore, Crocco’s calculations of cy, carried out only for 44, < 5 
and 7y,/7, < 2, can be supplemented by values of c; obtained 
with Sutherland’s viscosity law, even though different values of 
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cp and P may have been used. Values of (NP~'/#*R~'/2),, Caley. 
lated from Eq. (11) by using cy values obtained from references | 
and 20, are shown in Fig. 1. They are compared with the value 
of (NP-' t3R-'/2), given in references 1 and 20. The difference 
between the curves VD and C( YJ) or C( VD) shows the approxi. 
mate influence of the cp variation, while the difference between 
curves VJ and C( YJ) or C( VD) shows the influence of a decreas 
ing P. 

It is of further interest to note that Crocco’s values of cy for 
M, < 5 (Fig. 20, reference 18) slightly exceed the values of refer. 
ence 20, while the values based on a decreasing P! exceed thog 
of reference 18. Therefore, values of (NP~'/*R~'‘/*),, obtained 
with cy based on a slightly increasing P, would probably corre 
spond most closely with curves C(VD). Moreover, when Cp is 
considered variable, (NP~'/*R~'/?), appears to be fairly con- 
stant over an appreciable range of 44,; this corresponds to the 
usual engineering approximation. Also, at the larger values of 
M,, the effect of 7, diminishes. 
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Variation of local heat-transfer coefficient Parameter with Mach Number for laminar flow of constant pressure. 
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Effect of Flapping Hinge Offset on Rotor-Blade 
Stall 


J.P. Chawla 
Hughes Aircraft Company, 
1953 


} 


ilver Calif 


City 


February 24 


to some remarks on the subject by 

; Falabella, Jr., in the 
Readers’ Forum 1953, 
JOURNAI I should like only to point out the main differences 
between their approach and mine 


work by Meyer and Falabella are: 
a lifting rotor 


REPLY 
Meyer, IJr., and G 
the 


| HIS NOTE IS IN 
Messrs. John R 


Section of February, issue of the 


Significant features of the 


(1) The rotor is simply and is not inclined 
forward to provide a propulsive component of thrust as in a con 
ventional helicopter. A machine at an advance ratio of 0.60 with 
the tip-path plane of the rotor approximately horizontal is ob 
viously operating as a gyrodyne and not as a helicopter 

(2) The rotor is allowed to flap resulting in a hub moment that 
must be counterbalanced by a moment supplied by a control 
surface in order to trim the aircraft in flight 

I had taken the case of a conventional helicopter with normal 
values of Lock’s Inertia Number ¥ A value of y = 
by Meyer and Falabella is about a half to a third that of the cur 


It is recognized that with 


tas chosen 


rent practice even with tip-jet units 


hinge offset a reduction in the value of y is beneficial for rotor 


blade stall only if the blade is allowed to flap while in trimmed 
Consequently, the reduction in angle of attack at the tip 
1 of their 


flight 
note is due to a 


of the retreating blade shown in Fig 
the reduction in y; (2) the incorporation of 


the flapping of the rotor in trimmed 


combination of (1 
flapping hinge offset; (3 


flight; and (4) the operation of the machine as a gyrodyne with 

















he lieve 
stall 


the propulsive unit being other than the rotor itself. 1 
and (4) 


and (2); 


that the contribution of items (3) in reducing the 


preponderates over that of items (1 in other words, the 
blade mass constant and the flapping hinge offset play a secondary 


role 


ERRATA 
I should also like to submit corrections to my original note 


1) Inthe expression for \/,, Eq. (5a), the term 6;Cu should be 


preceded by a minus sign instead of a plus sign 


READERS’ FORUM 


Flow of Air in 


365 


: I I I 
(2 In Eqs. (5b), B = - e+ e 
4 2 } 
(3) In Eqs. (6a), EF = @&(A + Bu?) — 6Cu — AC 
, ‘ y | & l 
(4 In Eqs. (6b), Bo = ola (1 + wu?) — - (Op +A 
- > 


(5) In the numerator of Eq. (8), [¢/(1 + ¢)] should be pre 


ceded by a plus sign instead of a minus sign 


The effect of these errors is to decrease slightly the beneficial 
effect of flapping hinge offset on stall of the retreating blade tip 


shown in Fig. 2 of reference 1 
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Chawla, J. P., Effectof Flapping Hinge Offset on Rotor Blade Stall, Journal 


of the Aeronautical Sciences, Vol. 19, No. 7, p. 498, July, 1952 


Actuator Disc Theory for Compressible Flow 
and a Subsonic Correction for Propellers 


E. V. Laitone 
Associate Professor 
January 29, 1953 


_ ONE-DIMENSIONAL AXIAL-MOMENTUM THEORY can be ap 
plied to a propeller in subsonic flow if one assumes the con 
cept of the actuator dise creating a discontinuity in pressure and 
density. If we 
bility correction for small values of the free-stream Mach Number 
(Mo = Vo/Ae 


equations can easily be reduced to (e.g., using Eqs. (1) through (7) 


University of California at Berkeley 


are interested in only the first-order compressi 
and light disc loadings, then the axial-momentum 


of reference | 


= (1 1976 (1 
p 
Vi —, 
1 + ( l + VW ) 5; (2 
VW 2 
T/A T/A T/A 
5) I (1 Vy? | (3 
tg 7] | 2 po I 
ml = py VA p VoA |] + 6(1 VW? (4 
where 
Vv, = Vol + 6 (5 
| Vo (1 + 63) Vo (1 + 26; (6) 


using the same notation as in reference 1- namely, subscript (0 


indicates free-stream conditions an infinite distance ahead of 


propeller (1), and subscript (3) represents the conditions far 
behind in the propeller wake when p,; = p 

Now the simplest compressibility correction factor for the pro 
peller thrust is obtained by integrating Eq. (1.3) on page 251 of 
reference 2, assuming a constant 6;, constant chord, and constant 


angle of attack (a = 6 @) along the entire radius so 


2ra constant 


Cr cos } = (7) 
/ w \2 Vor%(1 + 6 frw\? 
v! , | 
ay \ a\* ay 
Cp sin @ 0 & 
Phen for light disc loadings (6; — 0), it is easily seen that 
Tcomp Li 3 
Tine 2 
1+ Me. . Vr V1 Mp? — Mit 
sin! 
Vy, V1 Mo? Up? 
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where 


Mo = Vo/ao, Mp = Rw/A 


and w is the actual rotational velocity of the propeller 

In order to apply Eq. (9) to finite disc loadings, it is easily shown 
that we must multiply it by Eq. (1) and replace J/) by 11, by 
means of Eq. (2) in order to get the higher order influence of the 
induced axial inflow (6, V»), considered uniform across the disc 
However, the question now arises whether this higher order term 
is justified in comparison with other neglected effects. For 
example, by using the same method as used in deriving Eq. (2.8) 
on page 194 of reference 2, it may be shown that for subsonic 
flow the first-order effect, due to the decreased enthalpy pro 
duced by the slipstream rotation, may be expressed as 


65 = (1/2)6; [1 + 61( MW? Mr?)| (10) 


Another term, also of the same order, is found to be that pro 
duced by the subsonic compressibility effect of the semi-infinite 
helical trailing vortex system accompanying the infinite number 
of propeller blades. This term may be evaluated to the first order 
by combining Eq. 48 of reference 3 with the procedure used in de 
riving Eq. (89) of reference 4 to show that 


l 
a — os l (11) 


for light disc loadings and Jo —~ 0. 

The important fact is now seen that Eqs. (1) and (11) have 
exactly the opposite influence on Eq. (9) as have Eqs. (2) and 
(10). Consequently, Eq. (9) is probably applicable to heavier 
disc loads in its present simplest form. The effect of varying in 
flow or varying chord can be shown to be a higher order effect 
that tends to cancel, leaving Eq. (9) as the ratio of the compres 
sibility effect for any propelier 

Eqs. (10) and (11) prove that the axial-momentum theory is 
not justified by itself alone. The question still arises, however, as 
to the limitations of the actuator disc concept in subsonic flow 
For example, it is easily shown that the so-called isentropic ‘‘ideal 
propeller” used in reference 1 is not a valid assumption for larger 
disc loadings. For Eq. (8) of reference 1 it is stated that the 
term due to p po exists because of the propeller profile drag and 
such additional losses. This, however, is not the case, since Eq. 
(8) of reference 1 corresponds to the required increase of entropy 
that must coexist with the discontinuous pressure rise assumed 
for the actuator disc. Therefore, one cannot assume both the 
Actually, Eqs 
(1) through (6), the first-order effect, should provide a better 


thin actuator disc and p; = py as in reference | 
approximation than reference 1. It is interesting to note that 
Eq. (4) approaches the correct limit for the mass-rate of flow for 


sonic conditions (J) —> 1 namely, 


m, =p VA = p Vy) rR? (12) 
The first-order effect of p py may be easily obtained from the 
energy equation and the first-order approximation for power re 
quired (7'V;) as 


showing that the effect of p; < po cannot be neglected unless 


|( po/ ps) 1] << (63)? (26, )? 
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Theory of Heat Transfer Through a Laminar 
Boundary Layer 


Myron Tribus 


Director of Icing Research, Engineering Research Institute 
University of Michigan, Ann Arbor, Mich. 
February 13, 1953 


i ee Readers’ Forum note by H. Schuh (JOURNAL OF THE 
AERONAUTICAL SCIENCES, Vol. 20, No. 2, p. 146, February, 
1953) appears to cover the same ground as the paper by M. J 
Lighthill on ‘‘Contributions to the Theory of Heat Transfer 
Through a Laminar Boundary Layer"’ (Proc. Roy. Soc. London, 
Series A, Vol. 202, pp. 359-377, 1950) 

The Lighthill results are somewhat more general than those 
given by Schuh 


Errata—‘‘Diffuser Efficiency of Free-Jet 
Supersonic Wind Tunnels at Variable Test 
Chamber Pressure’’+ 


Rudolf Hermann 


Professor, Department of Aeronautical Engineering, University of 
Minnesota, Minneapolis 


March 9, 1953 


iy THE REFERENCE ARTICLE,f the following errors were made 
(1) Eq. (29) should read 


2 M,* 


(2) On page 379, the fifth line following Eq. (42), ng must read 


nB- 
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Manuscripts: The original typewritten copy of the paper is 
desired, double or triple spaced on one side of white paper sheets, 
consecutively numbered. There should be wide margins to 
allow for the marking of directions to the printer. Correcting, 
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Sciences, 2 East 64th Street, New York 21, N.Y. All manu- 
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the paper is acceptable for publication. 


TitLes: The title of the paper should be brief. The name 
and initials of the author should be written as he prefers. The 
use of the full name of an author is advocated because of the fre- 
quent duplication of initials and surnames which sometimes 
makes it difficult to establish the identity of the author. This is 
particularly important for large annual indexing and abstracting 
services. All titles and degrees or honors are omitted. The 
name of the organization with which the author is associated 
should be placed after his name on a separate line. The date on 
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author’s title should be indicated in a footnote. 


SUMMARIES OR ApstrRaActTs: An abstract to be printed at the 
beginning should accompany each article. It should be adequate 
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major conclusions reached, since summaries in many cases con- 
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reference indexes. Abstracts printed in other journals, espe- 
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CAPTIONS AND LEGENDS: Legends or captions must accompany 
each drawing or photograph submitted. If written on the draw- 
ing or photograph, they should be placed below and well outside 
the part to be reproduced. Each table should have a caption 
such as Table 1, Table 2, Table 3, etc. Captions should be com- 
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